!----------------------------------------------------------------------!
! Optional advection scheme for CCATT-BRAMS/BRAMS models version 4.2+  ! 
! Based on Walcek, 2000 (JGR) and Walcek and Aleksic, 1998 (ATENV).    ! 
! The scheme is highly conservative, monotonic and keeps mass mixing   ! 
! ratio positive definite. 					       ! 
! Implemented by Saulo Freitas (saulo.freitas@cptec.inpe.br) @ Jun/2009! 
!----------------------------------------------------------------------!

MODULE monotonic_adv
  
  USE node_mod, only:        &
          	      ibcon, &  !INTENT(IN)
          	      mynum, &  !INTENT(IN)
          	      nodei0,   &  !INTENT(IN)
          	      nodej0,   &  !INTENT(IN)
          	      nodemyp,  &  !INTENT(IN)
          	      nodemxp,  &  !INTENT(IN)
          	      nodemzp   	!INTENT(IN)
  USE mem_grid, ONLY:        &
  	    	     dtlt,   & !INTENT(IN)\
		     time,   &
  	    	     ngrids, & !INTENT(IN)
  	    	     ngrid,  & !INTENT(IN)
	    	     dzt,    & !INTENT(IN)
	    	     dztn,   & !INTENT(IN)
            	     grid_g, & !INTENT(IN)	 
	    	     grid_g, & !INTENT(IN) 
	    	     naddsc, & !INTENT(IN) 
	             hw4       !INTENT(IN)	   
  
  USE mem_basic, ONLY: basic_g  !INTENT(IN) 
  
  USE micphys    ,  ONLY: level !INTENT(IN) 
  
  USE rconstants ,  ONLY: cp,p00,cv,rgas,cpi   !INTENT(IN)

  !lfr use mem_aer1, only : aerosol,num_scalar_aer_1st   !INTENT(IN)
  use mem_chem1, only : nspecies_transported        !INTENT(IN)

  use mem_scratch, only  : scratch  ! only scr1, inout
  
  use var_tables, only : scalar_tab & ! (var_p = IN, var_t = INOUT) 
                        ,num_scalar   ! (IN) 
			
  use advMessageMod, ONLY: SendMessageI,RecvMessageI,SendMessageJ,RecvMessageJ, &
                           newM2,newM3,newIa,newIz,newJa,newJz,nRecvI,nRecvJ,nSendI,nSendJ

  IMPLICIT NONE
  include "mpif.h"
  
  INTEGER:: advmnt

  INTEGER , PARAMETER :: ON=1,OFF=0
  
  INTEGER, PARAMETER :: GhostZoneLength=15

  INTEGER , PARAMETER :: use_true_density  = 0 ! 0= OFF, 1=ON
  
  !- for theoretical experiments
  INTEGER , PARAMETER :: theor_wind = 0        ! 0= OFF, 1=ON

  INTEGER           :: mnt_adv_jnitialized=0
  
  real, parameter :: c1 = cv/rgas, c2 = p00/rgas !c2 = p00*(cpi**c1)/rgas

  TYPE advmnt_vars 
     REAL,POINTER,DIMENSION  (:,:,:)  :: u3d 
     REAL,POINTER,DIMENSION  (:,:,:)  :: v3d 
     REAL,POINTER,DIMENSION  (:,:,:)  :: w3d 

     REAL,POINTER,DIMENSION  (:,:,:)  :: vc3d_jn 
     REAL,POINTER,DIMENSION  (:,:,:)  :: vc3d_out 
     REAL,POINTER,DIMENSION  (:,:,:)  :: vc3d_x 
     REAL,POINTER,DIMENSION  (:,:,:)  :: vc3d_y 

     REAL,POINTER,DIMENSION  (:,:,:)  :: dd0_3d 
     REAL,POINTER,DIMENSION  (:,:,:)  :: dd0_3du 
     REAL,POINTER,DIMENSION  (:,:,:)  :: dd0_3dv 
     REAL,POINTER,DIMENSION  (:,:,:)  :: dd0_3dw 

     REAL,POINTER,DIMENSION  (:,:,:)  :: den0_3d 
     REAL,POINTER,DIMENSION  (:,:,:)  :: den1_3d 
     REAL,POINTER,DIMENSION  (:,:,:)  :: den2_3d 
     REAL,POINTER,DIMENSION  (:,:,:)  :: den3_3d 

     REAL,POINTER,DIMENSION  (:,:)    :: dxtW 
     REAL,POINTER,DIMENSION  (:,:)    :: dytW 
     REAL,POINTER,DIMENSION  (:)      :: dztW 

  END TYPE advmnt_vars

  TYPE(advmnt_vars), ALLOCATABLE,DIMENSION(:)   :: advmnt_g,advmntNew_g
    
  PUBLIC :: advmnt_driver  ! Subroutine

  INTEGER, PARAMETER :: nvariables=2

  !nvariables=
  INTEGER, PARAMETER :: v_flux=1
  INTEGER, PARAMETER :: v_jmxmn=2
  
  
  INTEGER :: numberOfProcessors


  !type (watchtype) :: w(4)   ! This declares w to be a watch
  !type (watchtype) :: w1(6)   ! This declares w to be a watch
  
  INTEGER :: nSend_i,nSend_j,nRecv_i,nRecv_j
  INTEGER,PARAMETER :: bigdump=1
  REAL,ALLOCATABLE,DIMENSION(:,:,:) :: buffcomm
  INTEGER :: nRec_i,nSnd_i,nRec_j,nSnd_j
  INTEGER :: bufSendTotalLength_i,bufSendTotalLength_j
  INTEGER :: bufREcvTotalLength_i,bufRecvTotalLength_j

CONTAINS
  !----------------------------------------------------
             !call ADVECTc('V' ,mzp,mxp,myp,ia,iz,ja,jz,izu,jzv,mynum)
  SUBROUTINE advmnt_driver(varn,m1 ,m2 ,m3 ,ia,iz,ja,jz,izu,jzv,mynum)
  use node_mod, only : nodeibcon,nodeia,nodeiz,nodeja,nodejz,ipaths,iget_paths
    IMPLICIT NONE
    INTEGER , INTENT(IN) :: m1
    INTEGER , INTENT(IN) :: m2
    INTEGER , INTENT(IN) :: m3
    INTEGER , INTENT(IN) :: ia
    INTEGER , INTENT(IN) :: iz
    INTEGER , INTENT(IN) :: ja
    INTEGER , INTENT(IN) :: jz
    INTEGER , INTENT(IN) :: izu
    INTEGER , INTENT(IN) :: jzv
    INTEGER , INTENT(IN) :: mynum
    character(len=*),intent(IN) :: varn
    
    !--- local vars
    integer n,ng,mxyzp,i,j,procfile
    real, pointer :: scalarp, scalart
    
    character(4) :: ctime
        
    INTEGER :: iteste !lfr
    !PRINT *,'No driver da advect ',mynum;CALL flush(6)

    !call create_watch(w,name=(/"Waiting MPI "/))    ! Watches must be created before they are used
    !-scratch arrays initialization  and air density for using
    !- basic state densities
    
    !PRINT *, 'In advection........time=',time;CALL flush(6)
    
    iteste=0 !lfr
    
    IF(mnt_adv_jnitialized == OFF) THEN


       CALL initialize_advmnt(ngrids,nodemzp(mynum,:), &
                              nodemxp(mynum,:),nodemyp(mynum,:))


       do ng=1,ngrids
         CALL initialize_grid_spacings(ng,nodemzp(mynum,ng), &
	               nodemxp(mynum,ng),nodemyp(mynum,ng) &
                      ,grid_g(ng)%dxt	    &
                      ,grid_g(ng)%dyt	    &
 		      ,grid_g(ng)%fmapt     &
		      ,grid_g(ng)%rtgt      &
!
    		      ,advmnt_g(ng)%dxtW &
    		      ,advmnt_g(ng)%dytW &
		      ,advmnt_g(ng)%dztW )
       enddo
       
       if(use_true_density == OFF) then 
          do ng=1,ngrids
	    CALL initialize_densities(nodemzp(mynum,ng),&
	               nodemxp(mynum,ng),nodemyp(mynum,ng) &
                      , basic_g(ng)%dn0     &
		      , basic_g(ng)%dn0u    &
                      , basic_g(ng)%dn0v    &
                      ,advmnt_g(ng)%dd0_3d  &
    		      ,advmnt_g(ng)%dd0_3du &
    		      ,advmnt_g(ng)%dd0_3dv &
		      ,advmnt_g(ng)%dd0_3dw )
	   enddo
       endif
       
       mnt_adv_jnitialized= ON    
    END IF

    mxyzp=m1*m2*m3
    
    !- Advect  U, V, and W
    if (varn .eq. 'V' .or. varn .eq. 'ALL') then
      stop 'not using mnt to advect u,v,w'
    endif

    !- Advect  scalars
    
    !- get actual air densities, if using them instead of basic state fields
    if(use_true_density == ON) then
     call get_true_densities(m1,m2,m3,level &
                      , basic_g(ngrid)%rtp     &
		      , basic_g(ngrid)%rv      &
                      , basic_g(ngrid)%pp      &
		      , basic_g(ngrid)%pi0     &
                      , basic_g(ngrid)%theta   &
                      ,advmnt_g(ngrid)%dd0_3d  &
    		      ,advmnt_g(ngrid)%dd0_3du &
    		      ,advmnt_g(ngrid)%dd0_3dv &
		      ,advmnt_g(ngrid)%dd0_3dw )
    
    endif

 
    !- prepare wind velocities including map factors
    call prepare_winds(dtlt,m1,m2,m3,ia,iz,ja,jz     &
                      ,basic_g(ngrid)%uc  &
		      ,basic_g(ngrid)%up  &
                      ,basic_g(ngrid)%vc  &
                      ,basic_g(ngrid)%vp  &
		      ,basic_g(ngrid)%wc  &
                      ,basic_g(ngrid)%wp  &
!        
 		      ,grid_g(ngrid)%fmapui &
		      ,grid_g(ngrid)%fmapvi &
		      ,grid_g(ngrid)%rtgt   &
                      ,grid_g(ngrid)%rtgu   &
		      ,grid_g(ngrid)%rtgv   &
		      ,grid_g(ngrid)%f13t   &
                      ,grid_g(ngrid)%f23t   &
!
                      ,advmnt_g(ngrid)%u3d  &
    		      ,advmnt_g(ngrid)%v3d  &
    		      ,advmnt_g(ngrid)%w3d  )
                                             
 
    if(theor_wind == on) then
 
        call prepare_theor_winds(dtlt,m1,m2,m3,ia,iz,ja,jz     &
                      ,advmnt_g(ngrid)%u3d  &
    		      ,advmnt_g(ngrid)%v3d  &
    		      ,advmnt_g(ngrid)%w3d  &
                      ,grid_g(ngrid)%dxt    &
                      ,grid_g(ngrid)%dyt    &
                      ,advmnt_g(ngrid)%dd0_3d  &
    		      ,advmnt_g(ngrid)%dd0_3du &
    		      ,advmnt_g(ngrid)%dd0_3dv &
		      ,advmnt_g(ngrid)%dd0_3dw )
    endif
    !- prepare Walcek's air densities
    
    call get_Walceks_densities(dtlt,m1,m2,m3 &
                      ,advmnt_g(ngrid)%u3d  &
    		      ,advmnt_g(ngrid)%v3d  &
    		      ,advmnt_g(ngrid)%w3d  &
                      ,advmnt_g(ngrid)%dd0_3d  &
    		      ,advmnt_g(ngrid)%dd0_3du &
    		      ,advmnt_g(ngrid)%dd0_3dv &
		      ,advmnt_g(ngrid)%dd0_3dw &
                      ,advmnt_g(ngrid)%den0_3d &
    		      ,advmnt_g(ngrid)%den1_3d &
    		      ,advmnt_g(ngrid)%den2_3d &
		      ,advmnt_g(ngrid)%den3_3d &
    		      ,advmnt_g(ngrid)%dxtW &
    		      ,advmnt_g(ngrid)%dytW &
		      ,advmnt_g(ngrid)%dztW &
                      ,grid_g(ngrid)%dxt    &
                      ,grid_g(ngrid)%dyt    )

         !call advect_mnt_stalone()



    CALL copyInitFields(ngrid,m1,m2,m3)

!LFR>     CALL InitialFieldsUpdate(m1,m2,m3,newm2(ngrid),newm3(ngrid),ngrid,mynum, &
!LFR> 	   &
!LFR> 	   nRec_i, procRecv_i, tagRecv_i, &
!LFR> 	   iaRecv_i, izRecv_i, jaRecv_i, jzRecv_i, &
!LFR> 	   bufRecvStart_i, bufRecvLength_i, bufRecvTotalLength_i, &
!LFR>            &
!LFR> 	   nSnd_i, procSend_i, tagSend_i, &
!LFR> 	   iaSend_i, izSend_i, jaSend_i, sendMessageI(ngrid)%jz, &
!LFR> 	   bufSendStart_i, bufSendLength_i, bufSendTotalLength_i)
    
    CALL InitialFieldsUpdate(m1,m2,m3,newm2(ngrid),newm3(ngrid),ngrid,mynum, &
	   nrecvI(ngrid),RecvMessageI(ngrid)%proc,RecvMessageI(ngrid)%tag, &
	   RecvMessageI(ngrid)%ia,RecvMessageI(ngrid)%iz,RecvMessageI(ngrid)%ja,RecvMessageI(ngrid)%jz, &
	   RecvMessageI(ngrid)%start,RecvMessageI(ngrid)%mSize,TotalRecvI(ngrid), &
	   nSendI(ngrid),sendMessageI(ngrid)%proc,sendMessageI(ngrid)%tag, &
	   sendMessageI(ngrid)%ia,sendMessageI(ngrid)%iz,sendMessageI(ngrid)%ja,sendMessageI(ngrid)%jz, &
	   sendMessageI(ngrid)%start,sendMessageI(ngrid)%mSize,TotalSendI(ngrid))

!LFR>     CALL InitialFieldsUpdate(m1,m2,m3,newm2(ngrid),newm3(ngrid),ngrid,mynum, &
!LFR> 	   nRec_j, procRecv_j, tagRecv_j, &
!LFR> 	   iaRecv_j, izRecv_j, jaRecv_j, jzRecv_j, &
!LFR> 	   bufRecvStart_j, bufRecvLength_j, bufRecvTotalLength_j, &
!LFR> 	   nSnd_j, procSend_j, tagSend_j, &
!LFR> 	   iaSend_j, izSend_j, jaSend_j, jzSend_j, &
!LFR> 	   bufSendStart_j, bufSendLength_j, bufSendTotalLength_j)
	   
    CALL InitialFieldsUpdate(m1,m2,m3,newm2(ngrid),newm3(ngrid),ngrid,mynum, &
           nrecvJ(ngrid),RecvMessageJ(ngrid)%proc,RecvMessageJ(ngrid)%tag, &
	   RecvMessageJ(ngrid)%ia,RecvMessageJ(ngrid)%iz,RecvMessageJ(ngrid)%ja,RecvMessageJ(ngrid)%jz, &
	   RecvMessageJ(ngrid)%start,RecvMessageJ(ngrid)%mSize,TotalRecvJ(ngrid), &
           nSendJ(ngrid),sendMessageJ(ngrid)%proc,sendMessageJ(ngrid)%tag, &
	   sendMessageJ(ngrid)%ia,sendMessageJ(ngrid)%iz,sendMessageJ(ngrid)%ja,sendMessageJ(ngrid)%jz, &
	   sendMessageJ(ngrid)%start,sendMessageJ(ngrid)%mSize,TotalSendJ(ngrid))

STOP 188888880888888881

     !- ready to do advection, loop over all scalars
     do n=1,num_scalar(ngrid)     

!srf - somente para gases e aerossois
!     do n=num_scalar(ngrid) - NSPECIES_TRANSPORTED +1,num_scalar(ngrid)
!srf - somente para gases e aerossois
     !do n=num_scalar(ngrid) - NSPECIES_TRANSPORTED +1,num_scalar(ngrid)
        !print*,'nscalar1=',n
	!PRINT *,scalar_tab(n,ngrid)%name
	if (scalar_tab(n,ngrid)%name /= 'COP') cycle
        !print*,'nscalar2=',n
!srf - somente para gases e aerossois
!srf - somente para gases e aerossois
        !PRINT *,n; CALL flush(6)


 !srf-falta incluir controle para shaved_eta e aerossois com sedimentacao


      scalarp => scalar_tab(n,ngrid)%var_p
      scalart => scalar_tab(n,ngrid)%var_t

!for intel - descomente
      call atob(mxyzp,scalarp,advmnt_g(ngrid)%vc3d_jn) 
!LFR>       print *,'LFR  (.)(.) atob - passou ', mynum,n,ngrid,num_scalar(ngrid) ; call flush(6)

      call advect_mnt(ngrid,m1,m2,m3,ia,iz,ja,jz,dtlt &
!for pgi - descomente
!                     ,scalarp                 & ! initial mass mixing ratio (mmxr)
!for intel - descomente
 ,advmnt_g(ngrid)%vc3d_jn                     & ! initial mass mixing ratio (mmxr)
 ,advmnt_g(ngrid)%vc3d_out                    & ! final mmxr due advection
! ,advmntNew_g(ngrid)%vc3d_x			  &
! ,advmntNew_g(ngrid)%vc3d_y		     &
!-- 
		     ,advmntNew_g(ngrid)%u3d	 & ! 3D U wind
		    ,advmntNew_g(ngrid)%v3d	& ! 3D V wind
		    ,advmntNew_g(ngrid)%w3d	& ! 3D SIGMAZ wind
		     ,advmntNew_g(ngrid)%dd0_3d  & ! density 0
		     ,advmntNew_g(ngrid)%den0_3d & ! density D0; Eqs. 5  (JGR 2000)
		    ,advmntNew_g(ngrid)%den1_3d & ! density D1
		    ,advmntNew_g(ngrid)%den2_3d & ! density D2
		    ,advmntNew_g(ngrid)%den3_3d & ! density D3
		    ,advmntNew_g(ngrid)%dxtW	&
		    ,advmntNew_g(ngrid)%dytW	&
		    ,advmntNew_g(ngrid)%dztW	&

!		      ,advmnt_g(ngrid)%u3d     & ! 3D U wind
!		      ,advmnt_g(ngrid)%v3d     & ! 3D V wind
!		      ,advmnt_g(ngrid)%w3d     & ! 3D SIGMAZ wind
!		      ,advmnt_g(ngrid)%dd0_3d  & ! density 0
!		      ,advmnt_g(ngrid)%den0_3d & ! density D0; Eqs. 5  (JGR 2000)
!		      ,advmnt_g(ngrid)%den1_3d & ! density D1
!		      ,advmnt_g(ngrid)%den2_3d & ! density D2
!		      ,advmnt_g(ngrid)%den3_3d & ! density D3
!		      ,advmnt_g(ngrid)%dxtW    &
!		      ,advmnt_g(ngrid)%dytW    &
!		      ,advmnt_g(ngrid)%dztW    &

		     ,mynum,n)
      		     
      call advtndc(m1,m2,m3,ia,iz,ja,jz        &
!                  ,scalarp  ,scratch%scr1 (1)  &
                  ,scalarp  ,advmnt_g(ngrid)%vc3d_out  &
                  ,scalart  ,dtlt,mynum        )

        iteste=iteste+1
     enddo
!     PRINT * ,'Deallocate Comm',mynum;call flush(6)  
!     CALL CommDeAlloc()
!stop 333
 
  END SUBROUTINE advmnt_driver
 !----------------------------------------------------
  
  SUBROUTINE initialize_advmnt(ngrids, mmzp,mmxp,mmyp)
   implicit none 
   INTEGER , INTENT(IN) :: ngrids
   INTEGER , INTENT(IN) :: mmxp(ngrids)
   INTEGER , INTENT(IN) :: mmyp(ngrids)
   INTEGER , INTENT(IN) :: mmzp(ngrids)
   
   INTEGER :: ng
   
   !maxz=maxval(mmzp(1:ngrids))
   !maxx=maxval(mmxp(1:ngrids))
   !maxy=maxval(mmyp(1:ngrids))
   
   IF(ALLOCATED(advmnt_g)) THEN 
      PRINT *,'Error in initialize_advmnt, sub: radvc_mnt: advmnt_g already allocated!'
      RETURN
   END IF
   ALLOCATE (advmnt_g(ngrids))
   ALLOCATE (advmntNew_g(ngrids))

   do ng=1,ngrids
      ALLOCATE(advmnt_g(ng)%u3d    (mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%u3d=0.
      ALLOCATE(advmnt_g(ng)%v3d    (mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%v3d=0.
      ALLOCATE(advmnt_g(ng)%w3d    (mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%w3d=0.
      
      ALLOCATE(advmnt_g(ng)%dd0_3d (mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%dd0_3d =0.
      ALLOCATE(advmnt_g(ng)%dd0_3du(mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%dd0_3du=0.
      ALLOCATE(advmnt_g(ng)%dd0_3dv(mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%dd0_3dv=0.
      ALLOCATE(advmnt_g(ng)%dd0_3dw(mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%dd0_3dw=0.
      										      
      ALLOCATE(advmnt_g(ng)%den0_3d(mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%den0_3d=0.
      ALLOCATE(advmnt_g(ng)%den1_3d(mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%den1_3d=0.
      ALLOCATE(advmnt_g(ng)%den2_3d(mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%den2_3d=0.
      ALLOCATE(advmnt_g(ng)%den3_3d(mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%den3_3d=0.

      ALLOCATE(advmnt_g(ng)%dxtW(mmxp(ng),mmyp(ng))); advmnt_g(ng)%dxtW=0.
      ALLOCATE(advmnt_g(ng)%dytW(mmxp(ng),mmyp(ng))); advmnt_g(ng)%dytW=0.
      ALLOCATE(advmnt_g(ng)%dztW(mmzp(ng)))                  ; advmnt_g(ng)%dztW=0.


      ALLOCATE(advmnt_g(ng)%vc3d_jn  (mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%vc3d_jn =0.
      ALLOCATE(advmnt_g(ng)%vc3d_out (mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%vc3d_out=0.
      ALLOCATE(advmnt_g(ng)%vc3d_x   (mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%vc3d_x  =0.
      ALLOCATE(advmnt_g(ng)%vc3d_y   (mmzp(ng),mmxp(ng),mmyp(ng))); advmnt_g(ng)%vc3d_y  =0.

!-------

      ALLOCATE(advmntNew_g(ng)%u3d    (mmzp(ng),newM2(ng),newM3(ng))); advmntNew_g(ng)%u3d=0.
      ALLOCATE(advmntNew_g(ng)%v3d    (mmzp(ng),newM2(ng),newM3(ng))); advmntNew_g(ng)%v3d=0.
      ALLOCATE(advmntNew_g(ng)%w3d    (mmzp(ng),newM2(ng),newM3(ng))); advmntNew_g(ng)%w3d=0.
      
      ALLOCATE(advmntNew_g(ng)%dd0_3d (mmzp(ng),newM2(ng),newM3(ng))); advmntNew_g(ng)%dd0_3d =0.
      ALLOCATE(advmntNew_g(ng)%dd0_3du(mmzp(ng),newM2(ng),newM3(ng))); advmntNew_g(ng)%dd0_3du=0.
      ALLOCATE(advmntNew_g(ng)%dd0_3dv(mmzp(ng),newM2(ng),newM3(ng))); advmntNew_g(ng)%dd0_3dv=0.
      ALLOCATE(advmntNew_g(ng)%dd0_3dw(mmzp(ng),newM2(ng),newM3(ng))); advmntNew_g(ng)%dd0_3dw=0.
      										      
      ALLOCATE(advmntNew_g(ng)%den0_3d(mmzp(ng),newM2(ng),newM3(ng))); advmntNew_g(ng)%den0_3d=0.
      ALLOCATE(advmntNew_g(ng)%den1_3d(mmzp(ng),newM2(ng),newM3(ng))); advmntNew_g(ng)%den1_3d=0.
      ALLOCATE(advmntNew_g(ng)%den2_3d(mmzp(ng),newM2(ng),newM3(ng))); advmntNew_g(ng)%den2_3d=0.
      ALLOCATE(advmntNew_g(ng)%den3_3d(mmzp(ng),newM2(ng),newM3(ng))); advmntNew_g(ng)%den3_3d=0.

      ALLOCATE(advmntNew_g(ng)%dxtW(newM2(ng),newM3(ng))); advmntNew_g(ng)%dxtW=0.
      ALLOCATE(advmntNew_g(ng)%dytW(newM2(ng),newM3(ng))); advmntNew_g(ng)%dytW=0.
      ALLOCATE(advmntNew_g(ng)%dztW(mmzp(ng)))                  ; advmntNew_g(ng)%dztW=0.


      ALLOCATE(advmntNew_g(ng)%vc3d_jn  (mmzp(ng),newM2(ng),newM3(ng))); advmntNew_g(ng)%vc3d_jn =0.
      ALLOCATE(advmntNew_g(ng)%vc3d_out (mmzp(ng),newM2(ng),newM3(ng))); advmntNew_g(ng)%vc3d_out=0.
      ALLOCATE(advmntNew_g(ng)%vc3d_x   (mmzp(ng),newM2(ng),newM3(ng))); advmntNew_g(ng)%vc3d_x  =0.
      ALLOCATE(advmntNew_g(ng)%vc3d_y   (mmzp(ng),newM2(ng),newM3(ng))); advmntNew_g(ng)%vc3d_y  =0.

   enddo

   



  END SUBROUTINE initialize_advmnt
 !----------------------------------------------------
 
  SUBROUTINE initialize_densities(m1,m2,m3,dn0,dn0u,dn0v &
			    ,dd0_3d ,dd0_3du,dd0_3dv,dd0_3dw )
   IMPLICIT NONE 
   INTEGER , INTENT(IN) :: m1,m2,m3
   REAL,DIMENSION(m1,m2,m3),intent(IN) ::dn0,dn0u,dn0v
   REAL,DIMENSION(m1,m2,m3),intent(OUT)::dd0_3d,dd0_3du,dd0_3dv,dd0_3dw
   ! local var
   integer i,j,k
   
   dd0_3d (:,:,:)=  dn0 (:,:,:) 
   dd0_3du(:,:,:)=  dn0u(:,:,:)
   dd0_3dv(:,:,:)=  dn0v(:,:,:)
   do j = 1,m3
      do i = 1,m2
         do k = 1,m1-1
            dd0_3dw(k,i,j) = 0.5*(dn0(k,i,j) +dn0(k+1,i,j))
         enddo
	 dd0_3dw(m1,i,j)=dd0_3dw(m1-1,i,j)
   enddo;enddo
   

  END SUBROUTINE initialize_densities
 !----------------------------------------------------
 
  SUBROUTINE initialize_grid_spacings(ng,m1,m2,m3,dxt,dyt,fmapt,rtgt &
			    ,dxtW,dytW,dztW )
   IMPLICIT NONE 
   INTEGER , INTENT(IN) :: ng,m1,m2,m3
   REAL,DIMENSION(m2,m3),intent(IN) :: dxt,dyt,fmapt,rtgt
   
   REAL,DIMENSION(m2,m3),intent(OUT):: dxtW,dytW
   REAL,DIMENSION(m1),intent(OUT):: dztW
   ! local var
   integer i,j,k
   real rtgti
   
   do j = 1,m3
      do i = 1,m2
   	rtgti = 1. / rtgt(i,j)
        
	!- at init/rams_grid.f90:
        !     dxt(i,j)=fmapt(i,j)/(xmn(i,ngrid)-xmn(i-1,ngrid))  
        !     dyt(i,j)=fmapt(i,j)/(ymn(j,ngrid)-ymn(j-1,ngrid))

   	 dxtW(i,j) = 1./(dxt(i,j) * fmapt(i,j) * rtgti)
   	 dytW(i,j) = 1./(dyt(i,j) * fmapt(i,j) * rtgti)
   enddo;enddo
   do k = 1,m1
    !- at init/gridset.f90:	
    !  dztn(k,ifm) = 1. / (zmn(k,ifm) - zmn(k-1,ifm))
    ! Por que o Jacobiano nao depende de Z, o dztw depende somente 
    ! de z.
    !dztW(k,i,j) = 1./ ( dzt(k) * rtgti * fmapt(i,j)**2 )
     dztW(k)	 = 1./ ( dztn(k,ng) ) !
 
   enddo
  END SUBROUTINE initialize_grid_spacings
 !----------------------------------------------------
 
  SUBROUTINE get_true_densities(m1,m2,m3,level,rtp,rv,pp,pi0,theta &
			       ,dd0_3d ,dd0_3du,dd0_3dv,dd0_3dw )
   IMPLICIT NONE 
   INTEGER , INTENT(IN) :: m1,m2,m3,level
   REAL,DIMENSION(m1,m2,m3),intent(IN) ::rtp,rv,pp,pi0,theta
   REAL,DIMENSION(:,:,:),intent(OUT):: dd0_3d
   REAL,DIMENSION(m1,m2,m3),intent(OUT)  :: dd0_3du,dd0_3dv,dd0_3dw
   ! local var
   integer i,j,k
   real c3
   
   c3 = c2 * (cpi**c1)

   !- true air density at points "T"
   
   if( level == 0 ) then
     dd0_3d(:,:,:) = (c3/theta(:,:,:))*(pi0(:,:,:)+pp(:,:,:))**c1
   else
   do j = 1,m3
      do i = 1,m2
         do k = 1,m1
            dd0_3d(k,i,j) = (c3/theta(k,i,j))* (1. + rtp(k,i,j))/ &
	            (1. + 1.61*rv(k,i,j))*(pi0(k,i,j)+pp(k,i,j))**c1
         end do
      end do
   end do
   endif 

   !- true air density at points "U", "V" and "W":
   
   call fill_dn0uv(m1,m2,m3,dd0_3d,dd0_3du,dd0_3dv)

   do j = 1,m3
      do i = 1,m2
         do k = 1,m1-1
            dd0_3dw(k,i,j) = 0.5*(dd0_3d(k,i,j) + dd0_3d(k+1,i,j))
         enddo
	 dd0_3dw(m1,i,j)=dd0_3dw(m1-1,i,j)
   enddo;enddo

 END SUBROUTINE get_true_densities

 !----------------------------------------------------

 SUBROUTINE prepare_winds(dtlt,m1,m2,m3,ia,iz,ja,jz &
                            ,uc,up,vc,vp,wc,wp &
  		            ,fmapui &
		            ,fmapvi &
                            ,rtgt   &
                            ,rtgu   &
		            ,rtgv   &
		            ,f13t   &
                            ,f23t   &
			    ,u3d,v3d,w3d)
   
   implicit none 
   INTEGER , INTENT(IN) :: m1,m2,m3,ia,iz,ja,jz
   REAL, INTENT(IN) :: dtlt
   REAL,DIMENSION(m1,m2,m3),intent(INOUT) :: uc,up,vc,vp,wc,wp
!   REAL,DIMENSION(m1,m2,m3),intent(IN) :: uc,up,vc,vp,wc,wp ! TMP <<<<<<<<<<<<<<<<<<<<<<<
   REAL,DIMENSION(m2,m3)   ,intent(IN) :: rtgt,rtgu,rtgv,fmapui,fmapvi,f13t,f23t
   
   REAL,DIMENSION(m1,m2,m3),intent(OUT)::u3d,v3d,w3d
  
   !- local var
   real   dtlto2
   integer jm,jp,im,ip 
   integer i,j,k
   real :: c1,c2,rtgti

   dtlto2 = .5! * dtlt

  ! u3d, u3d, and w3d are input as the velocity components (averaged
  ! between past and current time levels) times dtlt.
   do j=1,m3
     do i = 1,m2
      do k = 1,m1
  
          w3d(k,i,j) = ( wc(k,i,j) + wp(k,i,j) )*dtlto2
          u3d(k,i,j) = ( uc(k,i,j) + up(k,i,j) )*dtlto2
          v3d(k,i,j) = ( vc(k,i,j) + vp(k,i,j) )*dtlto2
  
   enddo;enddo;enddo


  !return ! for pure cartesian coordinates

  ! here w3d is the cartesian vertical velocity 
  
  ! Add contribution to w3d from horiz winds crossing sloping sigma surfaces,
  ! and include 1/rtgt factor in w3d
  do j = 1,m3		 
     jm = max(1,j-1)
     jp = min(m3,j+1)
     do i = 1,m2
  	im = max(1,i-1)
  	ip = min(m2,i+1)
        rtgti = 1. / rtgt(i,j)
	
	do k = 1,m1-1
	    w3d(k,i,j) = ( (u3d(k,i,j) + u3d(k+1,i,j) + u3d(k,im,j) + u3d(k+1,im,j) ) * f13t(i,j)  &
	               +   (v3d(k,i,j) + v3d(k+1,i,j) + v3d(k,i,jm) + v3d(k+1,i,jm) ) * f23t(i,j)  &
		         ) * hw4(k)  &
	               + w3d(k,i,j) * rtgti

  	enddo
     enddo
  enddo
  ! here w3d is the sigma_z velocity 

  !- including map factors on U,V:
  do j = 1,m3
     do i = 1,m2
  	c1 = fmapui(i,j) * rtgu(i,j)
  	c2 = fmapvi(i,j) * rtgv(i,j)
  	do k = 1,m1-1
  	   u3d(k,i,j) = u3d(k,i,j) * c1 
  	   v3d(k,i,j) = v3d(k,i,j) * c2 
	   
	   
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
	   !u3d(k,i,j) = 0. !<<<<<<<<<<<<<<<<<<<<<<
	   
	   
	   !v3d(k,i,j) = -1. !<<<<<<<<<<<<<<<<<<<<<<

           !w3d(k,i,j) =	 0.
	   !wc(k,i,j)=w3d(k,i,j)
	   !uc(k,i,j)=u3d(k,i,j)
	   !vc(k,i,j)=v3d(k,i,j)
	   !wp(k,i,j)=w3d(k,i,j)
	   !up(k,i,j)=u3d(k,i,j)
	   !vp(k,i,j)=v3d(k,i,j)
	   
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  	enddo
     enddo
  enddo

  END SUBROUTINE prepare_winds
 !----------------------------------------------------

  SUBROUTINE get_Walceks_densities(dt,m1,m2,m3,u3d,v3d,w3d &
			    ,dd0_3d ,dd0_3du,dd0_3dv,dd0_3dw &
			    ,den0_3d,den1_3d,den2_3d,den3_3d &
			    ,dxtW,dytW,dztW,dxt,dyt)

   IMPLICIT NONE 
   !-in
   INTEGER , INTENT(IN) :: m1,m2,m3
   REAL    , INTENT(IN) :: dt
   REAL, DIMENSION(m1)      , INTENT(IN) :: dztW
   REAL, DIMENSION(m2,m3)   , INTENT(IN) :: dxtW,dytW

   REAL, DIMENSION(m2,m3)   , INTENT(IN) :: dxt,dyt

   REAL, DIMENSION(m1,m2,m3), INTENT(IN) :: u3d,v3d,w3d     &
			                   ,dd0_3du &
					   ,dd0_3dv,dd0_3dw
   REAL, DIMENSION(m1,m2,m3), INTENT(INOUT) :: dd0_3d 
   !-out
   REAL, DIMENSION(m1,m2,m3), INTENT(OUT):: den0_3d,den1_3d &
                                           ,den2_3d,den3_3d
   
   ! local var
   integer i,j,k
   
    DO  j=m3,2,-1
      DO  i=2,m2
        DO k = 2,m1
	

!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
!dd0_3d(k,i,j)=1.
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<


            den0_3d(k,i,j)=dd0_3d(k,i,j)



	    
!	    den1_3d(k,i,j)=den0_3d(k,i,j) - dt/dxt(i,j)*(d0         *u3d(k,i,j)- d0           *u3d(k,i-1,j ))
!	    den1_3d(k,i,j)=den0_3d(k,i,j) - dt*dxt(i,j)*(dd0_3du(k,i,j)*u3d(k,i,j)- dd0_3du(k,i-1,j)*u3d(k,i-1,j  ))
	    den1_3d(k,i,j)=den0_3d(k,i,j) - dt/dxtW(i,j)*(dd0_3du(k,i,j)*u3d(k,i,j)- dd0_3du(k,i-1,j)*u3d(k,i-1,j  ))

!           den2_3d(k,i,j)=den1_3d(k,i,j)- dt/dx*(d0*v3d(k,i,j)- d0*v3d(k,i,j-1 ))
!           den2_3d(k,i,j)=den1_3d(k,i,j)- dt*dyt(i,j)*(dd0_3dv(k,i,j)*v3d(k,i,j)- dd0_3dv(k,i,j-1)*v3d(k,i,j-1  ))
            den2_3d(k,i,j)=den1_3d(k,i,j)- dt/dytW(i,j)*(dd0_3dv(k,i,j)*v3d(k,i,j)- dd0_3dv(k,i,j-1)*v3d(k,i,j-1  ))
!
!           den3_3d(k,i,j)=den2_3d(k,i,j)-dt/dx         *(d0*w3d(k,i,j)-d0*w3d(  k-1,i,j))
!           den3_3d(k,i,j)=den2_3d(k,i,j)-dt*dzt(k)     *(dd0_3dw(k,i,j)*w3d(k,i,j)-dd0_3dw(k-1,i,j)*w3d(k-1,i,j))
            den3_3d(k,i,j)=den2_3d(k,i,j)-dt/dztW(k)*(dd0_3dw(k,i,j)*w3d(k,i,j)-dd0_3dw(k-1,i,j)*w3d(k-1,i,j))

       END DO !i
      END DO !j
    END DO !k  

 

 END SUBROUTINE get_Walceks_densities
 !----------------------------------------------------
 SUBROUTINE advect_mnt(ngrid,m1,m2,m3,ia,iz,ja,jz,dt,vc3d_in,vc3d_out &  
                      ,u3d,v3d,w3d,dd0,den0,den1,den2,den3 &
		      ,dxtW,dytW,dztW,mynum,n)
   
   use node_mod, only : nodei0,nodej0
   use mem_grid, only : time ! tmp
  IMPLICIT NONE
  INTEGER , INTENT(IN) :: m1,ngrid
  INTEGER , INTENT(IN) :: m2
  INTEGER , INTENT(IN) :: m3
  INTEGER , INTENT(IN) :: ia
  INTEGER , INTENT(IN) :: iz
  INTEGER , INTENT(IN) :: ja
  INTEGER , INTENT(IN) :: jz,n
  INTEGER , INTENT(IN) :: mynum
  REAL    , INTENT(IN) :: dt
  
  !REAL,DIMENSION(m1)      , INTENT(IN)  :: dzt
  !REAL,DIMENSION(m2,m3)   , INTENT(IN)  :: dxt,dyt,rtgt
  
  REAL,DIMENSION(m1),       INTENT(IN) :: dztW 
  REAL,DIMENSION(newM2(ngrid),newM3(ngrid))   , INTENT(IN) :: dxtW,dytW

  REAL,DIMENSION(m1,m2,m3), INTENT(IN) :: vc3d_in! initial mixing ratio 
  REAL,DIMENSION(m1,newM2(ngrid),newM3(ngrid)), INTENT(IN) :: u3d,v3d,w3d
  REAL,DIMENSION(m1,newM2(ngrid),newM3(ngrid)), intent(IN) :: dd0,den0,den1,den2,den3

  !REAL,DIMENSION(m1,m2,m3), INTENT(OUT) :: vc3d_x,vc3d_y
  REAL,DIMENSION(m1,m2,m3), INTENT(OUT) :: vc3d_out ! updated mixing ratio due advection

  !- local var
  REAL,DIMENSION(m1)               :: dxx  
  REAL,DIMENSION(m2,m3)            :: dxy  
  real masscon,initialmass,vol
  integer :: i,j,k,nrec,it
  
  CHARACTER(LEN=2) :: nf

  REAL,DIMENSION(m1,m2,m3) :: vc3dout_X,vc3dout_Y,vc3dout_Z
  REAL,DIMENSION(m1,newM2(ngrid),newM3(ngrid)) :: vc3X,vc3Y,vc3Z  

  INTEGER :: iBegin, iEnd, jBegin, jEnd
     
     iBegin=newIa(ngrid)-1
     iEnd=iBegin+m2-1
     jBegin=newJa(ngrid)-1
     jEnd=jBegin+m3-1
          

      !WRITE (71,*) m2,m3,newM2(ngrid),newM3(ngrid),newIa(ngrid),newJa(ngrid)
      !WRITE (71,*) iBegin,iEnd,jBegin,jEnd
      !CALL flush(71)

      vc3z(1:m1,iBegin:iEnd,jBegin:jEnd)=vc3d_in(1:m1,1:m2,1:m3)
!    vc3z=float(mynum)
    
! write(60+mynum,*) '1', vc3z(4,1:newM2(ngrid),7); call flush(60+mynum)
    
      call UpdateBorders(m1, newm2(ngrid), newm3(ngrid), vc3x, &
	 nrecvI(ngrid), RecvMessageI(ngrid)%proc, RecvMessageI(ngrid)%tag, &
	 RecvMessageI(ngrid)%ia, RecvMessageI(ngrid)%iz, RecvMessageI(ngrid)%ja, RecvMessageI(ngrid)%jz, &
	 RecvMessageI(ngrid)%start, RecvMessageI(ngrid)%mSize, TotalRecvI(ngrid), &
	 nSendI(ngrid), sendMessageI(ngrid)%proc, sendMessageI(ngrid)%tag, &
	 sendMessageI(ngrid)%ia, sendMessageI(ngrid)%iz, sendMessageI(ngrid)%ja, sendMessageI(ngrid)%jz, &
	 sendMessageI(ngrid)%start, sendMessageI(ngrid)%mSize, TotalSendI(ngrid))
!write(60+mynum,*) '2',vc3z(4,1:newM2(ngrid),7); call flush(60+mynum)

     call Advec3d_X(m1,newM2(ngrid),newM3(ngrid),2,newM2(ngrid)-1,2,newM3(ngrid)-1,vc3z,u3d,den0, &
                 den1,dt,dxtW,dd0,vc3x,mynum )
		 
!     call Advec3d_X(m1, newm2(ngrid), newm3(ngrid), &
!          vc3z, u3d, den0, den1, dt, dxtW, dd0, vc3x)
!srf tmp
!     vc3d_out(1:m1,1:m2,1:m3)=vc3x(1:m1,iBegin:iEnd,jBegin:jEnd)
!write(60+mynum,*) vc3x(4,1:newM2(ngrid),10); call flush(60+mynum)
 !    return
!srf tmp


     ! Do J-advection next over all I-k lines
     call UpdateBorders(m1, newm2(ngrid), newm3(ngrid), vc3z, &
          nrecvJ(ngrid), RecvMessageJ(ngrid)%proc, RecvMessageJ(ngrid)%tag, &
          RecvMessageJ(ngrid)%ia, RecvMessageJ(ngrid)%iz, RecvMessageJ(ngrid)%ja, RecvMessageJ(ngrid)%jz, &
          RecvMessageJ(ngrid)%start, RecvMessageJ(ngrid)%mSize, TotalRecvJ(ngrid), &
          nSendJ(ngrid), sendMessageJ(ngrid)%proc, sendMessageJ(ngrid)%tag, &
          sendMessageJ(ngrid)%ia, sendMessageJ(ngrid)%iz, sendMessageJ(ngrid)%ja, sendMessageJ(ngrid)%jz, &
          sendMessageJ(ngrid)%start, sendMessageJ(ngrid)%mSize, TotalSendJ(ngrid))

    call Advec3d_Y(m1,newM2(ngrid),newM3(ngrid),2,newM2(ngrid)-1,2,newM3(ngrid)-1,vc3x ,v3d,den1, &
                 den2,dt,dytW,dd0,vc3y )


!srf tmp
     vc3d_out(1:m1,1:m2,1:m3)=vc3y(1:m1,iBegin:iEnd,jBegin:jEnd)
     return
!srf tmp


!     call Advec3d_Y(m1, newm2(ngrid), newm3(ngrid), &
!          vc3x, v3d, den1, den2, dt, dytW, dd0, vc3y)

        !--- do k-advection next over all I-J lines
        ! vc3dout_Z = razao de mistura final (apos a adveccao nas 3 dimensoes)
             call Advec3d_Z(m1,newM2(ngrid),newM3(ngrid),2,newM2(ngrid)-1,2,newM3(ngrid)-1,vc3y, &
	                 w3d,den2,den3,dt,dztW,dd0,vc3z)
!LFR>         call Advec3d_Z(m1, newm2(ngrid), newm3(ngrid), &
!LFR>              vc3y, w3d, den2, den3, dt, dztW, dd0, vc3z)

   
     vc3d_out(1:m1,1:m2,1:m3)=vc3z(1:m1,iBegin:iEnd,jBegin:jEnd)

 END SUBROUTINE advect_mnt
 !----------------------------------------------------
 
 
  SUBROUTINE prepare_theor_winds(dtlt,m1,m2,m3,ia,iz,ja,jz &
                            ,u3d,v3d,w3d&
			    ,dxt,dyt &
			    ,dd0_3d ,dd0_3du,dd0_3dv,dd0_3dw )

   
   IMPLICIT NONE 
   INTEGER , INTENT(IN) :: m1,m2,m3,ia,iz,ja,jz
   REAL, INTENT(IN) :: dtlt
   REAL, DIMENSION(m2,m3)   , INTENT(IN) :: dxt,dyt

   REAL,DIMENSION(m1,m2,m3),intent(OUT)::u3d,v3d,w3d     &
                                        ,dd0_3d ,dd0_3du &
					,dd0_3dv,dd0_3dw 
  
   !- local var
   REAL   :: dtlto2
   INTEGER :: i,j,k
   REAL  :: ai0s  =  25.0
   REAL  :: aj0s  =  50.0
   REAL  :: umx   =   10.0
   REAL,PARAMETER :: pii   =   3.141592653589793
   REAL    :: umax  =   0.0
   REAL    :: anrev,curnt,rx,xa,ilop,iwndty,nrec
   
   dtlto2 =  10.!*dtlt
    
     !WRITE(6,*) ' Wind fields?  0-rotating; or  1-divergent winds'
     iwndty = 0  ! 0-rotating
     !iwndty = 1  ! 1-divergent winds

     IF(iwndty==1) ai0s= 50.5
     ilop= ai0s-21.  ! needed for printouts
     ! Define wind fields (rotation or divergent) and initial mixing ratios
     !  Cone at (25,50) for rotating winds; Cone at (50,50) divergent winds
     DO k = 1,m1
   	DO  j=m3,1,-1
   	   DO  i=1,m2
   	   
   	       dd0_3d (k,i,j)=1.
	       dd0_3du(k,i,j)=1. 
               dd0_3dv(k,i,j)=1.
	       dd0_3dw(k,i,j)=1. 

   	       u3d(k,i,j)= -2.*umx*(REAL(j)-REAL(110)/2.-.5)/REAL(110)*dtlto2
   	       v3d(k,i,j)=  2.*umx*(REAL(i)-REAL(100)/2.-.5)/REAL(100)*dtlto2
 	       w3d(k,i,j)= 0.                                        *dtlto2
   	   
   	      IF(iwndty==1) THEN
   		 xa=pii/25.
   		 IF(J>0) u3d(k,i,j)=umx*SIN(xa*REAL(i))*SIN(xa*(REAL(j)))
   		 IF(I>0) v3d(k,i,j)=umx*COS(xa*(REAL(i)-.5))*COS(XA*(REAL(j)+.5))
   	      END IF
   	   
   	   
   	      umax= MAX(ABS(u3d(k,i,j)),ABS(v3d(k,i,j)),umax)
   	      rx= SQRT((REAL(i)-ai0s)**2.+(REAL(j)-aj0s)**2.)

   	   END DO !i
   	END DO !j
     END DO !k  

  END SUBROUTINE prepare_theor_winds
  
   subroutine UpdateBorders(m1, m2, m3, field, &
	nRecv, procRecv, tagRecv, iaRecv, izRecv, jaRecv, jzRecv, &
	bufRecvStart, bufRecvLength, bufRecvTotalLength, &
	nSend, procSend, tagSend, iaSend, izSend, jaSend, jzSend, &
	bufSendStart, bufSendLength, bufSendTotalLength)
     integer, intent(in) :: m1
     integer, intent(in) :: m2
     integer, intent(in) :: m3
     real,    intent(inout) :: field(m1,m2,m3)
     integer, intent(in) :: nRecv
     integer, intent(in) :: procRecv(nRecv)
     integer, intent(in) :: tagRecv(nRecv)
     integer, intent(in) :: iaRecv(nRecv)
     integer, intent(in) :: izRecv(nRecv)
     integer, intent(in) :: jaRecv(nRecv)
     integer, intent(in) :: jzRecv(nRecv)
     integer, intent(in) :: bufRecvStart(nRecv)
     integer, intent(in) :: bufRecvLength(nRecv)
     integer, intent(in) :: bufRecvTotalLength
     integer, intent(in) :: nSend
     integer, intent(in) :: procSend(nSend)
     integer, intent(in) :: tagSend(nSend)
     integer, intent(in) :: iaSend(nSend)
     integer, intent(in) :: izSend(nSend)
     integer, intent(in) :: jaSend(nSend)
     integer, intent(in) :: jzSend(nSend)
     integer, intent(in) :: bufSendStart(nSend)
     integer, intent(in) :: bufSendLength(nSend)
     integer, intent(in) :: bufSendTotalLength
 
 
     integer :: i, j, iCnt, i1, i2, cnt, iRecv, iSend, ierr
     real, allocatable :: bufRecv(:)
     integer :: reqRecv(nRecv)
     real, allocatable :: bufSend(:)
     integer :: reqSend(nSend)
     integer :: status(MPI_STATUS_SIZE)
     integer :: status2(MPI_STATUS_SIZE,nSend)
 
     allocate(bufRecv(bufRecvTotalLength))
     do iRecv = 1, nRecv
	call MPI_Irecv(bufRecv(bufRecvStart(iRecv)), bufRecvLength(iRecv), &
	     MPI_REAL, &
	     procRecv(iRecv), tagRecv(iRecv), MPI_COMM_WORLD, &
	     reqRecv(iRecv), ierr)
     end do
 
     allocate(bufSend(bufSendTotalLength))
     do iSend = 1, nSend
	i1 = bufSendStart(iSend)
	iCnt = bufSendStart(iSend)
	i2 = bufSendLength(iSend)
	do j = jaSend(iSend), jzSend(iSend)
	   do i = iaSend(iSend), izSend(iSend)
	      bufSend(iCnt:iCnt+m1-1) = field(1:m1,i,j)
	      iCnt = iCnt+m1
	   end do
	end do
	call MPI_Isend(bufSend(i1), i2, MPI_REAL, &
	     procSend(iSend), tagSend(iSend), MPI_COMM_WORLD, &
	     reqSend(iSend), ierr)
     end do
 
     do cnt = 1, nRecv
	call MPI_Waitany(nRecv, reqRecv, iRecv, status, ierr)
	i1 = bufRecvStart(iRecv)
	iCnt = bufRecvStart(iRecv)
	i2 = bufRecvLength(iRecv)
	do j = jaRecv(iRecv), jzRecv(iRecv)
	   do i = iaRecv(iRecv), izRecv(iRecv)
	      field(1:m1,i,j) = bufRecv(iCnt:iCnt+m1-1)
	      iCnt = iCnt+m1
	   end do
	end do
     end do
 
     call MPI_Waitall(nSend, reqSend, status2, ierr)
 
     deallocate(bufSend)
     deallocate(bufRecv)
   end subroutine UpdateBorders
  
  SUBROUTINE InitialFieldsUpdate(m1,m2,m3,Nm2,Nm3,ng,mynum, &
	nRec, procRecv, tagRecv, iaRecv, izRecv, jaRecv, jzRecv, &
	bufRecvStart, bufRecvLength, bufRecvTotalLength, &
	nSnd, procSend, tagSend, iaSend, izSend, jaSend, jzSend, &
	bufSendStart, bufSendLength, bufSendTotalLength)
  
     integer, intent(in) :: m1
     integer, intent(in) :: m2,nm2
     integer, intent(in) :: m3,nm3,ng,mynum
     integer, intent(in) :: nRec
     integer, intent(in) :: procRecv(nRec)
     integer, intent(in) :: tagRecv(nRec)
     integer, intent(in) :: iaRecv(nRec)
     integer, intent(in) :: izRecv(nRec)
     integer, intent(in) :: jaRecv(nRec)
     integer, intent(in) :: jzRecv(nRec)
     integer, intent(in) :: bufRecvStart(nRec)
     integer, intent(in) :: bufRecvLength(nRec)
     integer, intent(in) :: bufRecvTotalLength
     integer, intent(in) :: nSnd
     integer, intent(in) :: procSend(nSnd)
     integer, intent(in) :: tagSend(nSnd)
     integer, intent(in) :: iaSend(nSnd)
     integer, intent(in) :: izSend(nSnd)
     integer, intent(in) :: jaSend(nSnd)
     integer, intent(in) :: jzSend(nSnd)
     integer, intent(in) :: bufSendStart(nSnd)
     integer, intent(in) :: bufSendLength(nSnd)
     integer, intent(in) :: bufSendTotalLength

     REAL,DIMENSION (m1,nm2,nm3) :: l_dxtW,l_dytW
     INTEGER :: i,j,k
     
    IF(bufSendTotalLength==0 .or. bufRecvTotalLength==0) RETURN
    
     DO i=1,nm2
     	DO j=1,nm3
     	   DO k=1,m1
     	      l_dxtW(k,i,j)=advmntNew_g(ng)%dxtW(i,j)
     	      l_dytW(k,i,j)=advmntNew_g(ng)%dytW(i,j)
     	   END DO
     	END DO
     END DO


    call UpdateBorders(m1, nm2, nm3,advmntNew_g(ng)%u3d, &
     nRec, procRecv, tagRecv, &
     iaRecv, izRecv, jaRecv, jzRecv, &
     bufRecvStart, bufRecvLength, bufRecvTotalLength, &
     nSnd, procSend, tagSend, &
     iaSend, izSend, jaSend, jzSend, &
     bufSendStart, bufSendLength, bufSendTotalLength)
!    WRITE (*,FMT='("Em init. Sou o ",I2.2,", na pos. ", I2.2)') myNum,1; CALL flush(6)

    call UpdateBorders(m1, nm2, nm3,advmntNew_g(ng)%v3d, &
     nRec, procRecv, tagRecv, &
     iaRecv, izRecv, jaRecv, jzRecv, &
     bufRecvStart, bufRecvLength, bufRecvTotalLength, &
     nSnd, procSend, tagSend, &
     iaSend, izSend, jaSend, jzSend, &
     bufSendStart, bufSendLength, bufSendTotalLength)
!    WRITE (*,FMT='("Em init. Sou o ",I2.2,", na pos. ", I2.2)') myNum,2; CALL flush(6)
     
    call UpdateBorders(m1, nm2, nm3,advmntNew_g(ng)%w3d, &
     nRec, procRecv, tagRecv, &
     iaRecv, izRecv, jaRecv, jzRecv, &
     bufRecvStart, bufRecvLength, bufRecvTotalLength, &
     nSnd, procSend, tagSend, &
     iaSend, izSend, jaSend, jzSend, &
     bufSendStart, bufSendLength, bufSendTotalLength)
!    WRITE (*,FMT='("Em init. Sou o ",I2.2,", na pos. ", I2.2)') myNum,3; CALL flush(6)
    call UpdateBorders(m1, nm2, nm3,advmntNew_g(ng)%vc3d_jn, &
     nRec, procRecv, tagRecv, &
     iaRecv, izRecv, jaRecv, jzRecv, &
     bufRecvStart, bufRecvLength, bufRecvTotalLength, &
     nSnd, procSend, tagSend, &
     iaSend, izSend, jaSend, jzSend, &
     bufSendStart, bufSendLength, bufSendTotalLength)
!    WRITE (*,FMT='("Em init. Sou o ",I2.2,", na pos. ", I2.2)') myNum,4; CALL flush(6)

    call UpdateBorders(m1, nm2, nm3,advmntNew_g(ng)%vc3d_out, &
     nRec, procRecv, tagRecv, &
     iaRecv, izRecv, jaRecv, jzRecv, &
     bufRecvStart, bufRecvLength, bufRecvTotalLength, &
     nSnd, procSend, tagSend, &
     iaSend, izSend, jaSend, jzSend, &
     bufSendStart, bufSendLength, bufSendTotalLength)
!    WRITE (*,FMT='("Em init. Sou o ",I2.2,", na pos. ", I2.2)') myNum,5; CALL flush(6)

    call UpdateBorders(m1, nm2, nm3,advmntNew_g(ng)%vc3d_x, &
     nRec, procRecv, tagRecv, &
     iaRecv, izRecv, jaRecv, jzRecv, &
     bufRecvStart, bufRecvLength, bufRecvTotalLength, &
     nSnd, procSend, tagSend, &
     iaSend, izSend, jaSend, jzSend, &
     bufSendStart, bufSendLength, bufSendTotalLength)
!    WRITE (*,FMT='("Em init. Sou o ",I2.2,", na pos. ", I2.2)') myNum,6; CALL flush(6)

    call UpdateBorders(m1, nm2, nm3,advmntNew_g(ng)%vc3d_y, &
     nRec, procRecv, tagRecv, &
     iaRecv, izRecv, jaRecv, jzRecv, &
     bufRecvStart, bufRecvLength, bufRecvTotalLength, &
     nSnd, procSend, tagSend, &
     iaSend, izSend, jaSend, jzSend, &
     bufSendStart, bufSendLength, bufSendTotalLength)
!    WRITE (*,FMT='("Em init. Sou o ",I2.2,", na pos. ", I2.2)') myNum,7; CALL flush(6)

    call UpdateBorders(m1, nm2, nm3,advmntNew_g(ng)%dd0_3d, &
     nRec, procRecv, tagRecv, &
     iaRecv, izRecv, jaRecv, jzRecv, &
     bufRecvStart, bufRecvLength, bufRecvTotalLength, &
     nSnd, procSend, tagSend, &
     iaSend, izSend, jaSend, jzSend, &
     bufSendStart, bufSendLength, bufSendTotalLength)
!    WRITE (*,FMT='("Em init. Sou o ",I2.2,", na pos. ", I2.2)') myNum,8; CALL flush(6)

    call UpdateBorders(m1, nm2, nm3,advmntNew_g(ng)%dd0_3du, &
     nRec, procRecv, tagRecv, &
     iaRecv, izRecv, jaRecv, jzRecv, &
     bufRecvStart, bufRecvLength, bufRecvTotalLength, &
     nSnd, procSend, tagSend, &
     iaSend, izSend, jaSend, jzSend, &
     bufSendStart, bufSendLength, bufSendTotalLength)
!    WRITE (*,FMT='("Em init. Sou o ",I2.2,", na pos. ", I2.2)') myNum,9; CALL flush(6)

    call UpdateBorders(m1, nm2, nm3,advmntNew_g(ng)%dd0_3dv, &
     nRec, procRecv, tagRecv, &
     iaRecv, izRecv, jaRecv, jzRecv, &
     bufRecvStart, bufRecvLength, bufRecvTotalLength, &
     nSnd, procSend, tagSend, &
     iaSend, izSend, jaSend, jzSend, &
     bufSendStart, bufSendLength, bufSendTotalLength)
!    WRITE (*,FMT='("Em init. Sou o ",I2.2,", na pos. ", I2.2)') myNum,10; CALL flush(6)

    call UpdateBorders(m1, nm2, nm3,advmntNew_g(ng)%dd0_3dw, &
     nRec, procRecv, tagRecv, &
     iaRecv, izRecv, jaRecv, jzRecv, &
     bufRecvStart, bufRecvLength, bufRecvTotalLength, &
     nSnd, procSend, tagSend, &
     iaSend, izSend, jaSend, jzSend, &
     bufSendStart, bufSendLength, bufSendTotalLength)
!    WRITE (*,FMT='("Em init. Sou o ",I2.2,", na pos. ", I2.2)') myNum,11; CALL flush(6)

    call UpdateBorders(m1, nm2, nm3,advmntNew_g(ng)%den0_3d, &
     nRec, procRecv, tagRecv, &
     iaRecv, izRecv, jaRecv, jzRecv, &
     bufRecvStart, bufRecvLength, bufRecvTotalLength, &
     nSnd, procSend, tagSend, &
     iaSend, izSend, jaSend, jzSend, &
     bufSendStart, bufSendLength, bufSendTotalLength)
!    WRITE (*,FMT='("Em init. Sou o ",I2.2,", na pos. ", I2.2)') myNum,12; CALL flush(6)

    call UpdateBorders(m1, nm2, nm3,advmntNew_g(ng)%den1_3d, &
     nRec, procRecv, tagRecv, &
     iaRecv, izRecv, jaRecv, jzRecv, &
     bufRecvStart, bufRecvLength, bufRecvTotalLength, &
     nSnd, procSend, tagSend, &
     iaSend, izSend, jaSend, jzSend, &
     bufSendStart, bufSendLength, bufSendTotalLength)
!    WRITE (*,FMT='("Em init. Sou o ",I2.2,", na pos. ", I2.2)') myNum,13; CALL flush(6)

    call UpdateBorders(m1, nm2, nm3,advmntNew_g(ng)%den2_3d, &
     nRec, procRecv, tagRecv, &
     iaRecv, izRecv, jaRecv, jzRecv, &
     bufRecvStart, bufRecvLength, bufRecvTotalLength, &
     nSnd, procSend, tagSend, &
     iaSend, izSend, jaSend, jzSend, &
     bufSendStart, bufSendLength, bufSendTotalLength)
!    WRITE (*,FMT='("Em init. Sou o ",I2.2,", na pos. ", I2.2)') myNum,14; CALL flush(6)

    call UpdateBorders(m1, nm2, nm3,advmntNew_g(ng)%den3_3d, &
     nRec, procRecv, tagRecv, &
     iaRecv, izRecv, jaRecv, jzRecv, &
     bufRecvStart, bufRecvLength, bufRecvTotalLength, &
     nSnd, procSend, tagSend, &
     iaSend, izSend, jaSend, jzSend, &
     bufSendStart, bufSendLength, bufSendTotalLength)
!    WRITE (*,FMT='("Em init. Sou o ",I2.2,", na pos. ", I2.2)') myNum,15; CALL flush(6)

    call UpdateBorders(m1, nm2, nm3,l_dxtW, &
     nRec, procRecv, tagRecv, &
     iaRecv, izRecv, jaRecv, jzRecv, &
     bufRecvStart, bufRecvLength, bufRecvTotalLength, &
     nSnd, procSend, tagSend, &
     iaSend, izSend, jaSend, jzSend, &
     bufSendStart, bufSendLength, bufSendTotalLength)
!    WRITE (*,FMT='("Em init. Sou o ",I2.2,", na pos. ", I2.2)') myNum,16; CALL flush(6)

    call UpdateBorders(m1, nm2, nm3,l_dytW, &
     nRec, procRecv, tagRecv, &
     iaRecv, izRecv, jaRecv, jzRecv, &
     bufRecvStart, bufRecvLength, bufRecvTotalLength, &
     nSnd, procSend, tagSend, &
     iaSend, izSend, jaSend, jzSend, &
     bufSendStart, bufSendLength, bufSendTotalLength)
!    WRITE (*,FMT='("Em init. Sou o ",I2.2,", na pos. ", I2.2)') myNum,17; CALL flush(6)

      DO i=1,nm2
     	DO j=1,nm3
     	      advmntNew_g(ng)%dxtW(i,j)=l_dxtW(1,i,j)
     	      advmntNew_g(ng)%dytW(i,j)=l_dytW(1,i,j)
     	END DO
     END DO

  END SUBROUTINE InitialFieldsUpdate

  SUBROUTINE CopyInitFields(ngrid,m1,m2,m3)
     INTEGER, INTENT(IN) :: ngrid,m1,m2,m3
     INTEGER :: iBegin,iEnd,jBegin,jEnd
     
     iBegin=newIa(ngrid)-1
     iEnd=newIz(ngrid)+1
     jBegin=newJa(ngrid)-1
     jEnd=newJz(ngrid)+1
     
!     IF(bigDump==1) 
!   WRITE (*,fmt='(A,5(I3.3,1X))') 'Copy limits: ',mynum,iBegin,iEnd,jBegin,jEnd
     
   advmntNew_g(ngrid)%u3d(1:m1,iBegin:iEnd,jBegin:jEnd)=advmnt_g(ngrid)%     u3d(1:m1,1:m2,1:m3) 
   advmntNew_g(ngrid)%v3d(1:m1,iBegin:iEnd,jBegin:jEnd)=advmnt_g(ngrid)%     v3d(1:m1,1:m2,1:m3) 
   advmntNew_g(ngrid)%w3d(1:m1,iBegin:iEnd,jBegin:jEnd)=advmnt_g(ngrid)%     w3d(1:m1,1:m2,1:m3) 

   !advmntNew_g(ngrid)%vc3d_jn(1:m1,iBegin:iEnd,jBegin:jEnd)=advmnt_g(ngrid)% vc3d_jn(1:m1,1:m2,1:m3)
   !advmntNew_g(ngrid)%vc3d_out(1:m1,iBegin:iEnd,jBegin:jEnd)=advmnt_g(ngrid)%vc3d_out(1:m1,1:m2,1:m3)
   advmntNew_g(ngrid)%vc3d_x(1:m1,iBegin:iEnd,jBegin:jEnd)=advmnt_g(ngrid)%  vc3d_x(1:m1,1:m2,1:m3) 
   advmntNew_g(ngrid)%vc3d_y(1:m1,iBegin:iEnd,jBegin:jEnd)=advmnt_g(ngrid)%  vc3d_y(1:m1,1:m2,1:m3) 

   advmntNew_g(ngrid)%dd0_3d(1:m1,iBegin:iEnd,jBegin:jEnd)=advmnt_g(ngrid)%  dd0_3d(1:m1,1:m2,1:m3) 
   advmntNew_g(ngrid)%dd0_3du(1:m1,iBegin:iEnd,jBegin:jEnd)=advmnt_g(ngrid)% dd0_3du(1:m1,1:m2,1:m3)
   advmntNew_g(ngrid)%dd0_3dv(1:m1,iBegin:iEnd,jBegin:jEnd)=advmnt_g(ngrid)% dd0_3dv(1:m1,1:m2,1:m3)
   advmntNew_g(ngrid)%dd0_3dw(1:m1,iBegin:iEnd,jBegin:jEnd)=advmnt_g(ngrid)% dd0_3dw(1:m1,1:m2,1:m3)

   advmntNew_g(ngrid)%den0_3d(1:m1,iBegin:iEnd,jBegin:jEnd)=advmnt_g(ngrid)% den0_3d(1:m1,1:m2,1:m3)
   advmntNew_g(ngrid)%den1_3d(1:m1,iBegin:iEnd,jBegin:jEnd)=advmnt_g(ngrid)% den1_3d(1:m1,1:m2,1:m3)
   advmntNew_g(ngrid)%den2_3d(1:m1,iBegin:iEnd,jBegin:jEnd)=advmnt_g(ngrid)% den2_3d(1:m1,1:m2,1:m3)
   advmntNew_g(ngrid)%den3_3d(1:m1,iBegin:iEnd,jBegin:jEnd)=advmnt_g(ngrid)% den3_3d(1:m1,1:m2,1:m3)

   advmntNew_g(ngrid)%dxtW(iBegin:iEnd,jBegin:jEnd)=advmnt_g(ngrid)%    dxtW(1:m2,1:m3) 
   advmntNew_g(ngrid)%dytW(iBegin:iEnd,jBegin:jEnd)=advmnt_g(ngrid)%    dytW(1:m2,1:m3) 
   advmntNew_g(ngrid)%dztW(1:m1)=advmnt_g(ngrid)%    dztW(1:m1)   

  END SUBROUTINE CopyInitFields

 
SUBROUTINE Advec3d_X(m1,m2,m3, ia,iz,ja,jz,q0,u,den0,den1,dt,dxx,dd0,qn,mynum)
!-------------------------
  ! This subroutine calculates change in mixing ratio (Q0) during time
  !  step DT due to advection along a grid IDIM in length. Mixing ratios
  !  from host code (C) are loaded into Q0 array, which is updated to QN.
  !  Velocities (U) and fluxes (FLUX) are specified at cell FACES, having
  !  dimensions 0:IDIM. U, Q0, QN, DXX and FLUX indices defined here:
  !  Densities at beg, end time (DEN0, DEN1) defined in HOST CODE  
  !
  ! I grid->   |  1  |  2  |  I-1  |   I  |..   ..|  IDIM  | <- host grid
  ! U-array-> u(0)  u(1)  u(2)   u(i-1)  u(i)           u(IDIM)
  ! C-array->  | C(1)| C(2)| C(I-1)| C(I) |..   ..| C(IDIM)| mixing ratio
  ! DXX-arry-> | Dx1 | Dx2 | DxI-1 | DxI  |..   ..| DxIDIM |
  ! Density->  | Dd1 | Dd2 | DdI-1 | DdI  |..   ..| DdIDIM |
  !                 Q0 defined along 0 - IDIM+1 cells:
  !    |       | QN  | QN  |  QN   |  QN  |       |   QN   |        |
  
  !    |   Q0--|-Q0--|-Q0--|--Q0 --|--Q0--|..   ..|-- Q0 --|--Q0    |
  !    |    0  | 1   |  2  | I-1   |  I   |       |  IDIM  | IDIM+1 |
  !   lower BC |             <---   Q0 grid   --->         | upper BC
  !           Boundary conditions are stored in Q0 cells 0 & IDIM+1 
  !
  !  Input to this subroutine, provided in common /sub/, and the calling
  !  arguments to this subroutine:
  !     IDIM - #of grid cells being updated
  !     Q0(0:IDIM+1)- Initial mixing ratio along 1-D array, with two
  !                 additional boundary value mixing ratios padded into the
  !                 0th and IDIM+1 cell locations
  !     U(0:IDIM)- velocities BETWEEN grid cells (at the "higher-I" edges of
  !                each grid cell in the array, units consistent with DX, DT
  !     DEN0(IDIM)- Initial fluid density, which needs to be updated during
  !                 multi-dimensional calculations, as noted in Calling code
  !     DEN1(IDIM)- Updated fluid density, which needs to be updated during
  !                 multi-dimensional calculations, as noted in calling code
  !     DT-         time step- units consistent with U
  !     DXX(IDIM)-  Grid cell length along advection direction, Units
  !                   consistent with DT and U
  !     DD0(0:IDIM)- Initial fluid density flowing BETWEEN each grid cell
  !                  (remains constant for all dimensions at the initial
  !                  fluid density of the 1st dimension of a 2-3 D calculation
  !               one can use UPSTREAM density here (DD0(I)= RHO0(I) if u>0
  !               or DD0(I)= RHO0(I+1) if u<0) where RHO0 is the initial
  !               fluid density at the beginning of the 1st dimensional
  !               advection step of a 2 or 3 D advection calculation done one
  !               step at a time 
  !
  !  Output of this subroutine is an updated mixing ratio array QN(IDIM)
  !

  use node_mod, only : nodeibcon,nodeia,nodeiz,nodeja,nodejz, &
                       ipaths,iget_paths,nodei0,nodej0,nmachs
       
  use mem_grid, only:  ngrids, nnxp,nnyp

  IMPLICIT none  
 
  INCLUDE "i8.h"
 
  INTEGER,INTENT(IN)                        :: m1,m2,m3, ia,iz,ja,jz,mynum
  REAL   ,INTENT(inout),DIMENSION(m1,m2,m3) :: q0
  REAL   ,INTENT(inout),DIMENSION(m1,m2,m3)    :: u
  REAL   ,INTENT(inout),DIMENSION(m1,m2,m3)    :: den0
  REAL   ,INTENT(inout),DIMENSION(m1,m2,m3)    :: den1
  REAL   ,INTENT(in)                        :: dt
  REAL   ,INTENT(in),DIMENSION(m2,m3)       :: dxx
  REAL   ,INTENT(inout),DIMENSION(m1,m2,m3)    :: dd0
  REAL   ,INTENT(OUT),DIMENSION(m1,m2,m3)   :: qn

  REAL,DIMENSION(m1,m2,m3)    :: flux
  REAL,DIMENSION(m1,m2,m3)    :: vcmax
  REAL,DIMENSION(m1,m2,m3)    :: vcmin
  LOGICAL,DIMENSION(m1,m2,m3) :: imxmn
  REAL,PARAMETER                       :: zr0=0.0
  
  INTEGER :: idime
  INTEGER :: i,j,k
  REAL    :: cf,cf1,ck1,ck2,x1,x1n

  integer, parameter :: debug =0
  integer, parameter :: on = 1
  
  INTEGER :: ii,ji,ii0,ji0,ie,je,ie0,je0,ipos,iia,iiz,nvar
  INTEGER,DIMENSION(nmachs*nmachs*ngrids*2*nvariables) :: isend_req,irecv_req
  
  
  integer :: nf
  
  !nvariables=
  !             1- flux
  REAL,DIMENSION(nvariables,m1) :: varSend
  REAL,DIMENSION(nvariables,m1) :: varRecv  
  
  nf=80+mynum
 idime = m2 ! x dir

  qn=q0

! WRITE (nf,FMT='("Idime= ",I3.3)') idime
! WRITE (nf,FMT='("m2=",I3.3," ia=",I3.3," iz=",I3.3," Map(m2)=",I3.3," Map(ia)=",I3.3," Map(iz)=",I3.3)') &
!       m2,ia,iz,m2+nodei0(mynum,ngrid),ia+nodei0(mynum,ngrid),iz+nodei0(mynum,ngrid);CALL flush(nf)
 
  imxmn=.false.
  ! imxmn(idime+1)=.false.

 DO j=ja,jz

    DO k=2,m1-1 
       ! Update mixing ratios and limit Fluxes going UP where u>0
       !  First assume upstream flux at edge of domain
       IF(u(k,1,j)>=zr0) flux(k,1,j)= q0(k,1,j)*u(k,1,j)*dt*dd0(k,1,j)
    END DO
  END DO 
 
 DO j=ja,jz
 
  ! Identify local max and min, specify mixing ratio limits at new time
  DO  i=2,idime-1 ! ia,iz-1 or 1,iz-1

     DO k=2,m1-1 
				  
        imxmn(k,i,j)=q0(k,i,j)>=max(q0(k,i-1,j),q0(k,i+1,j)) .or. & !=true if local
              q0(k,i,j)<=min(q0(k,i-1,j),q0(k,i+1,j))       !       extrema
        ck1= q0(k,i,j)
        ck2= q0(k,i,j)
        if(u(k,i,j  )< zr0) ck1= q0(k,i+1,j)
        if(u(k,i-1,j)>=zr0) ck2= q0(k,i-1,j)
        vcmax(k,i,j)= max( q0(k,i,j), ck1, ck2 )                      ! Eq-7
        vcmin(k,i,j)= min( q0(k,i,j), ck1, ck2 )                      ! Eq-7
        !  VCMAX and VCMIN are the absolute physical limits to the
        !     mixing ratio at t+dt. If these limits are ever violated,
        !     non-monotonic (oscillatory) behavior in solution results
        !

        IF(u(k,i,j)<zr0) CYCLE
     
        IF(u(k,i-1,j)<zr0) THEN
           flux(k,i,j)= q0(k,i,j)*u(k,i,j)*dt*dd0(k,i,j)    !  outflow-only cell     
        ELSE                              !      use upstream
           x1= dt*u(k,i,j)/dxx(i,j)               ! Courant number
           x1n= (1.-x1)*(q0(k,i+1,j)-q0(k,i-1,j))/4.
           !
           ! First, estimate mixing ratio in outflowing fluid (Cf)
           cf= q0(k,i,j) + x1n                                       !Eq-4a
           !
           !   Check to see if there is a peak (min) upwind and/or
           !    downwind of cell face
           IF(imxmn(k,i-1,j)) cf= q0(k,i,j) +MAX(1.5,1.2  +.6 *x1)*x1n   !Eq-10b
           IF(imxmn(k,i+1,j)) cf= q0(k,i,j) +       (1.75 -.45*x1)*x1n   !Eq-10a
           !        CF= Q0(k,i,j) + 5.*X1N   ! uncomment this line for "full sharp"
           !
           !   Limit Cf to be between mixing ratio on either side of edge
           !      where flux is being calculated
           cf1= MIN( MAX( cf, MIN(q0(k,i,j),q0(k,i+1,j))  ), MAX(q0(k,i,j),q0(k,i+1,j)) )
           !- for debug purposes only
	   if(debug==on)        CF1= CF     ! This statement IGNORES monotonic limitations.
           !                       you should uncomment this line and run your
           !                       advection calculation wth constant initial
           !                       mixing ratios everywhere. If you have properly
           !                       implemented this subroutine constant mixing
           !                       ratios should be maintained
           ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
        
	   !
           !   Calculate mixing ratio at new time, but limit to physically
           !    reasonable values
           qn(k,i,j)= MAX(vcmin(k,i,j),MIN(vcmax(k,i,j),          &   !eq-3&8
                    (q0(k,i,j)*den0(k,i,j)-x1*cf1*dd0(k,i,j)+flux(k,i-1,j)/dxx(i,j))/den1(k,i,j) ))

	   ! for debug purposes only
           if(debug==on)  QN(k,i,j)=(Q0(k,i,j)*DEN0(k,i,j)-X1*CF1*DD0(k,i,j)+FLUX(k,I-1,j)/DXX(I,J))/DEN1(k,i,j)
           !            This statement IGNORES monotonic limitations.
           !                       you should uncomment this line and run your
           !                       advection calculation wth constant initial
           !                       mixing ratios everywhere. If you have properly
	   !                       implemented this subroutine constant mixing
           !                       ratios should be maintiained
           ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
           !   print*,'i q= ', i,qn(i),line,q0(i),q0(i-1),q0(i+1)
	   !       if(qn(i) < 99.999) pause

	   !
           !   Re-calculate OUTFLOWING flux before moving on to next cell
           !    Flux = CF1*X1*DD0 but it must be adjusted if a monotonic limit
           !    is encountered.
           flux(k,i,j)= dxx(i,j)*(q0(k,i,j)*den0(k,i,j) - qn(k,i,j)*den1(k,i,j)) + flux(k,i-1,j)
        END IF                                                  !Eq-9a
     END DO
     
  END DO
END DO
  !
  ! If periodic boundary conditions are assumed, it is necessary
  !   to recalculate the updated mixing ratio at cell 1 if there
  !   is inflow to that cell from the boundary between IDIM and 1
  !   Here these statements are commented out, but should be uncommented
  !   if this subroutine is needed for periodic boundary conditions,
  !   and then one of the calling arguements to the subroutine is IPERIOD
  !   which is set to "1" if you assume period boundary conditions
  !      IF(IPERIOD==1) THEN
  !        IF(U(IDIM-1)>=ZR0.AND.U(IDIM)>=ZR0)
  !     &  QN(1)=(Q0(1)*DEN0(1)-FLUX(1)/DXX(1)+FLUX(IDIM)/DXX(1))/DEN1(1)
  !      END IF 
  !
  ! Update mixing ratios and limit Fluxes going DOWN where u<0
  !  The logic of this loop through the grid line is identical
  !  to the "DO 10" Loop above, only you start at the highest I
  !  edge and work backwards to I=1
  !

 DO j=ja,jz

    DO k=2,m1-1       
       IF(u(k,idime-1,j)<zr0) flux(k,idime-1,j)= &
    	       q0(k,idime,j)*u(k,idime-1,j)*dt*dd0(k,idime-1,j)
    END DO
  END DO 


  DO j=ja,jz
   
   !LFR-> from iz to ia
   DO i=idime-1,2,-1 !iz,ia,-1
     ! print*,'part 2'
     DO k=2,m1-1 
        
        IF(u(k,i-1,j)>=zr0) THEN           ! Inflow-only cell
        
	   IF(u(k,i,j)<zr0) qn(k,i,j)=  MAX(  vcmin(k,i,j),   MIN(   vcmax(k,i,j),&
                  (q0(k,i,j)*den0(k,i,j)-flux(k,i,j)/dxx(i,j) + flux(k,i-1,j)/dxx(i,j))/den1(k,i,j) ))
           ELSE
              x1=  dt*ABS(u(k,i-1,j))/dxx(i,j)     ! Courant number
              x1n= (1.-x1)*(q0(k,i-1,j)-q0(k,i+1,j))/4.
              cf= q0(k,i,j) + x1n                                       !Eq-4b
              IF(imxmn(k,i+1,j)) cf= q0(k,i,j) +MAX(1.5,1.2  +.6 *x1)*x1n   !Eq-10b
              IF(imxmn(k,i-1,j)) cf= q0(k,i,j) +   (1.75 -.45*x1)*x1n       !Eq-10a
              cf1= MIN( MAX( cf, MIN(q0(k,i,j),q0(k,i-1,j)) ), MAX(q0(k,i,j),q0(k,i-1,j)) )
              !- for debug purposes only
	      if(debug==on)  CF1= CF	! This statement IGNORES monotonic limitations.
              !			   you should uncomment this line and run your
              !			   advection calculation wth constant initial
              !			   mixing ratios everywhere. If you have properly
              !			   implemented this subroutine constant mixing
              !			   ratios should be maintained
              ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
              
	      IF(u(k,i,j)>=zr0) cf1= q0(k,i,j)     ! outflow-only cell upstream
              qn(k,i,j)= MAX(  vcmin(k,i,j),  MIN(   vcmax(k,i,j), 	  &   !Eq-3&8
                    (q0(k,i,j)*den0(k,i,j)-flux(k,i,j)/dxx(i,j)-x1*cf1*dd0(k,i-1,j))/den1(k,i,j) ))
              !print*,'i q= ', i,qn(i),line,q0(i),q0(i-1),q0(i+1)
          
	      !- for debug purposes only
              if(debug==on) QN(k,i,j)=(Q0(k,i,j)*DEN0(k,i,j)-FLUX(k,I,j)/DXX(I,J)-X1*CF1*DD0(k,I-1,j))/DEN1(k,i,j)
              !		This statement IGNORES monotonic limitations.
              !			   you should uncomment this line and run your
              !			   advection calculation wth constant initial
              !			   mixing ratios everywhere. If you have properly
              !			   implemented this subroutine constant mixing
	      !			   ratios should be maintiained
              ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
              !  print*,'i q= ', i,qn(i),line,q0(i),q0(i-1),q0(i+1)
              !  print*,FLUX(I),x1,cf1!,dd0(i-1),den1(i)
              !  print*,q0(i)-FLUX(I)-x1*cf1!,dd0(i-1),den1(i)
	      !      if(qn(i) < 99.999) pause

	      flux(k,i-1,j)=dxx(i,j)*(qn(k,i,j)*den1(k,i,j) - q0(k,i,j)*den0(k,i,j)) + flux(k,i,j)!Eq-9b
        END IF                                                  
     END DO

  END DO

  !END DO
  !
  ! If periodic boundary conditions are assumed, it is necessary
  !   to recalculate the updated mixing ratio at cell IDIM if there
  !   is inflow to that cell from the boundary between IDIM and 1
  !   Here these statements are commented out, but should be uncommented
  !   if this subroutine is needed for periodic boundary conditions,
  !   and then one of the calling arguements to the subroutine is IPERIOD
  !   which is set to "1" if you assume period boundary conditions
  !      IF(IPERIOD==1) THEN
  !      IF(U(1).LT.ZR0.AND.U(IDIM).LT.ZR0)
  !xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
  !     &  QN(IDIM)=(Q0(IDIM)*DEN0(IDIM)-FLUX(0)/DXX(IDIM)+FLUX(IDIM-1)/
  !     &                DXX(IDIM))/DEN1(IDIM)
  !      END IF
  !
ENDDO !- big loop y-z


END SUBROUTINE Advec3d_X

!---------------------------------------------------------------------------

SUBROUTINE Advec3d_Y(m1,m2,m3, ia,iz,ja,jz,q0,u,den0,den1,dt,dxx,dd0,qn)
!-------------------------
  ! This subroutine calculates change in mixing ratio (Q0) during time
  !  step DT due to advection along a grid IDIM in length. Mixing ratios
  !  from host code (C) are loaded into Q0 array, which is updated to QN.
  !  Velocities (U) and fluxes (FLUX) are specified at cell FACES, having
  !  dimensions 0:IDIM. U, Q0, QN, DXX and FLUX indices defined here:
  !  Densities at beg, end time (DEN0, DEN1) defined in HOST CODE  
  !
  ! I grid->   |  1  |  2  |  I-1  |   I  |..   ..|  IDIM  | <- host grid
  ! U-array-> u(0)  u(1)  u(2)   u(i-1)  u(i)           u(IDIM)
  ! C-array->  | C(1)| C(2)| C(I-1)| C(I) |..   ..| C(IDIM)| mixing ratio
  ! DXX-arry-> | Dx1 | Dx2 | DxI-1 | DxI  |..   ..| DxIDIM |
  ! Density->  | Dd1 | Dd2 | DdI-1 | DdI  |..   ..| DdIDIM |
  !                 Q0 defined along 0 - IDIM+1 cells:
  !    |       | QN  | QN  |  QN   |  QN  |       |   QN   |        |
  !    |   Q0--|-Q0--|-Q0--|--Q0 --|--Q0--|..   ..|-- Q0 --|--Q0    |
  !    |    0  | 1   |  2  | I-1   |  I   |       |  IDIM  | IDIM+1 |
  !   lower BC |             <---   Q0 grid   --->         | upper BC
  !           Boundary conditions are stored in Q0 cells 0 & IDIM+1 
  !
  !  Input to this subroutine, provided in common /sub/, and the calling
  !  arguments to this subroutine:
  !     IDIM - #of grid cells being updated
  !     Q0(0:IDIM+1)- Initial mixing ratio along 1-D array, with two
  !                 additional boundary value mixing ratios padded into the
  !                 0th and IDIM+1 cell locations
  !     U(0:IDIM)- velocities BETWEEN grid cells (at the "higher-I" edges of
  !                each grid cell in the array, units consistent with DX, DT
  !     DEN0(IDIM)- Initial fluid density, which needs to be updated during
  !                 multi-dimensional calculations, as noted in Calling code
  !     DEN1(IDIM)- Updated fluid density, which needs to be updated during
  !                 multi-dimensional calculations, as noted in calling code
  !     DT-         time step- units consistent with U
  !     DXX(IDIM)-  Grid cell length along advection direction, Units
  !                   consistent with DT and U
  !     DD0(0:IDIM)- Initial fluid density flowing BETWEEN each grid cell
  !                  (remains constant for all dimensions at the initial
  !                  fluid density of the 1st dimension of a 2-3 D calculation
  !               one can use UPSTREAM density here (DD0(I)= RHO0(I) if u>0
  !               or DD0(I)= RHO0(I+1) if u<0) where RHO0 is the initial
  !               fluid density at the beginning of the 1st dimensional
  !               advection step of a 2 or 3 D advection calculation done one
  !               step at a time 
  !
  !  Output of this subroutine is an updated mixing ratio array QN(IDIM)
  !
  use node_mod, only : nodeibcon,nodeia,nodeiz,nodeja,nodejz, &
                       ipaths,iget_paths,nodei0,nodej0,nmachs
       
  use mem_grid, only:  ngrids, nnxp,nnyp


  IMPLICIT none  
 
  INCLUDE "i8.h"
  
  INTEGER,INTENT(IN)                   :: m1,m2,m3, ia,iz,ja,jz
  INTEGER                  :: idime
  REAL   ,INTENT(inout),DIMENSION(m1,m2,m3) :: q0
  REAL   ,INTENT(inout),DIMENSION(m1,m2,m3)   :: u
  REAL   ,INTENT(inout),DIMENSION(m1,m2,m3)     :: den0
  REAL   ,INTENT(inout),DIMENSION(m1,m2,m3)     :: den1
  REAL   ,INTENT(in)                   :: dt
  REAL   ,INTENT(in),DIMENSION(m2,m3)     :: dxx
  REAL   ,INTENT(inout),DIMENSION(m1,m2,m3)   :: dd0
  REAL   ,INTENT(OUT),DIMENSION(m1,m2,m3)    :: qn

  REAL,DIMENSION(m1,m2,m3)                 :: flux
  REAL,DIMENSION(m1,m2,m3)                   :: vcmax
  REAL,DIMENSION(m1,m2,m3)                   :: vcmin
  LOGICAL,DIMENSION(m1,m2,m3)            :: imxmn
  REAL,PARAMETER                       :: zr0=0.0
  
  INTEGER :: i,j,k
  REAL    :: cf,cf1,ck1,ck2,x1,x1n

  integer, parameter :: debug =0
  integer, parameter :: on = 1
  !
  INTEGER :: ii,ji,ii0,ji0,ie,je,ie0,je0,ipos,iia,iiz,nvar
  INTEGER,DIMENSION(nmachs*nmachs*ngrids*2*nvariables) :: isend_req,irecv_req
  
  
  integer :: nf
  
  !nvariables=
  !             1- q
  !             2- flux
  REAL,DIMENSION(nvariables,m1) :: varSend
  REAL,DIMENSION(nvariables,m1) :: varRecv  
  
  nf=80+mynum

  idime = m3 ! y dir

  qn= q0
  imxmn=.false.
  ! imxmn(idime+1)=.false.

 DO i=ia,iz

    ! Update mixing ratios and limit Fluxes going UP where u>0
    !  First assume upstream flux at edge of domain
    DO k=2,m1-1   
       IF(u(k,i,1)>=zr0) flux(k,i,1)= q0(k,i,1)*u(k,i,1)*dt*dd0(k,i,1)
    END DO
      
  END DO 

 !- big loop y-z
 DO i=ia,iz
   
  ! Identify local max and min, specify mixing ratio limits at new time
  DO  j=2,idime-1 ! ja,jz
     DO k=2,m1-1   
        imxmn(k,i,j)=q0(k,i,j)>=max(q0(k,i,j-1),q0(k,i,j+1)) .or. & !=true if local
                     q0(k,i,j)<=min(q0(k,i,j-1),q0(k,i,j+1))	    !	    extrema
        ck1= q0(k,i,j)
        ck2= q0(k,i,j)
        if(u(k,i,j  )< zr0) ck1= q0(k,i,j+1)
        if(u(k,i,j-1)>=zr0) ck2= q0(k,i,j-1)
        vcmax(k,i,j)= max( q0(k,i,j), ck1, ck2 )                      ! Eq-7
        vcmin(k,i,j)= min( q0(k,i,j), ck1, ck2 )                      ! Eq-7
        !  VCMAX and VCMIN are the absolute physical limits to the
        !	mixing ratio at t+dt. If these limits are ever violated,
        !	non-monotonic (oscillatory) behavior in solution results
        !print*,'part 1'
        IF(u(k,i  ,j)<zr0) CYCLE
     
        IF(u(k,i,j-1)<zr0) THEN
           flux(k,i,j)= q0(k,i,j)*u(k,i,j)*dt*dd0(k,i,j)    !  outflow-only cell
     
        ELSE                              !      use upstream
           x1= dt*u(k,i,j)/dxx(i,j)               ! Courant number
           x1n= (1.-x1)*(q0(k,i,j+1)-q0(k,i,j-1))/4.
           !
           ! First, estimate mixing ratio in outflowing fluid (Cf)
           cf= q0(k,i,j) + x1n                                       !Eq-4a
           !
           !   Check to see if there is a peak (min) upwind and/or
           !    downwind of cell face
           IF(imxmn(k,i,j-1)) cf= q0(k,i,j) +MAX(1.5,1.2  +.6 *x1)*x1n   !Eq-10b
           IF(imxmn(k,i,j+1)) cf= q0(k,i,j) +       (1.75 -.45*x1)*x1n   !Eq-10a
           !        CF= Q0(k,i,j) + 5.*X1N   ! uncomment this line for "full sharp"
           !
           !   Limit Cf to be between mixing ratio on either side of edge
           !      where flux is being calculated
           cf1= MIN( MAX( cf, MIN(q0(k,i,j),q0(k,i,j+1))  ), MAX(q0(k,i,j),q0(k,i,j+1)) )
           !- for debug purposes only
	   if(debug==on)        CF1= CF     ! This statement IGNORES monotonic limitations.
           !                       you should uncomment this line and run your
           !                       advection calculation wth constant initial
           !                       mixing ratios everywhere. If you have properly
           !                       implemented this subroutine constant mixing
           !                       ratios should be maintained
           ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
        
	   !
           !   Calculate mixing ratio at new time, but limit to physically
           !    reasonable values
           qn(k,i,j)= MAX(  vcmin(k,i,j),   MIN(   vcmax(k,i,j),          &   !eq-3&8
                 (q0(k,i,j)*den0(k,i,j)-x1*cf1*dd0(k,i,j)+flux(k,i,j-1)/dxx(i,j))/den1(k,i,j) ))

	   ! for debug purposes only
           if(debug==on)  QN(k,i,j)=(Q0(k,i,j)*DEN0(k,i,j)-X1*CF1*DD0(k,i,j)+FLUX(k,i,J-1)/DXX(I,J))/DEN1(k,i,j)
           !            This statement IGNORES monotonic limitations.
           !                       you should uncomment this line and run your
           !                       advection calculation wth constant initial
           !                       mixing ratios everywhere. If you have properly
	   !                       implemented this subroutine constant mixing
           !                       ratios should be maintiained
           ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
           !   print*,'i q= ', i,qn(i),line,q0(i),q0(i-1),q0(i+1)
	   !       if(qn(i) < 99.999) pause

	   !
           !   Re-calculate OUTFLOWING flux before moving on to next cell
           !    Flux = CF1*X1*DD0 but it must be adjusted if a monotonic limit
           !    is encountered.
           flux(k,i,j)= dxx(i,j)*(q0(k,i,j)*den0(k,i,j) - qn(k,i,j)*den1(k,i,j)) + flux(k,i,j-1)
        END IF                                                  !Eq-9a
     END DO
  END DO
  !
  ! If periodic boundary conditions are assumed, it is necessary
  !   to recalculate the updated mixing ratio at cell 1 if there
  !   is inflow to that cell from the boundary between IDIM and 1
  !   Here these statements are commented out, but should be uncommented
  !   if this subroutine is needed for periodic boundary conditions,
  !   and then one of the calling arguements to the subroutine is IPERIOD
  !   which is set to "1" if you assume period boundary conditions
  !      IF(IPERIOD==1) THEN
  !        IF(U(IDIM-1)>=ZR0.AND.U(IDIM)>=ZR0)
  !     &  QN(1)=(Q0(1)*DEN0(1)-FLUX(1)/DXX(1)+FLUX(IDIM)/DXX(1))/DEN1(1)
  !      END IF 
  !
  ! Update mixing ratios and limit Fluxes going DOWN where u<0
  !  The logic of this loop through the grid line is identical
  !  to the "DO 10" Loop above, only you start at the highest I
  !  edge and work backwards to I=1
  !
  END DO
  
 DO i=ia,iz

    DO k=2,m1-1   
       IF(u(k,i,idime-1)<zr0) flux(k,i,idime-1)=q0(k,i,idime)*u(k,i,idime-1)*dt*dd0(k,i,idime-1)
    END DO
  END DO 

  DO i=ia,iz
  
  DO j=idime-1,2,-1 !jz,ja,-1
     DO k=2,m1-1     
        IF(u(k,i,j-1)>=zr0) THEN           ! Inflow-only cell        
	   IF(u(k,i,j)<zr0) qn(k,i,j)=  MAX(  vcmin(k,i,j),   MIN(   vcmax(k,i,j),&
                  (q0(k,i,j)*den0(k,i,j)-flux(k,i,j)/dxx(i,j) + flux(k,i,j-1)/dxx(i,j))/den1(k,i,j) ))
           ELSE
              x1=  dt*ABS(u(k,i,j-1))/dxx(i,j)     ! Courant number
              x1n= (1.-x1)*(q0(k,i,j-1)-q0(k,i,j+1))/4.
              cf= q0(k,i,j) + x1n                                       !Eq-4b
              IF(imxmn(k,i,j+1)) cf= q0(k,i,j) +MAX(1.5,1.2  +.6 *x1)*x1n   !Eq-10b
              IF(imxmn(k,i,j-1)) cf= q0(k,i,j) +   (1.75 -.45*x1)*x1n       !Eq-10a
              cf1= MIN( MAX( cf, MIN(q0(k,i,j),q0(k,i,j-1)) ), MAX(q0(k,i,j),q0(k,i,j-1)) )
              !- for debug purposes only
	      if(debug==on)  CF1= CF	! This statement IGNORES monotonic limitations.
              !			   you should uncomment this line and run your
              !			   advection calculation wth constant initial
              !			   mixing ratios everywhere. If you have properly
              !			   implemented this subroutine constant mixing
              !			   ratios should be maintained
              ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
              
	      IF(u(k,i,j)>=zr0) cf1= q0(k,i,j)     ! outflow-only cell upstream
              qn(k,i,j)= MAX(  vcmin(k,i,j),  MIN(   vcmax(k,i,j), 	  &   !Eq-3&8
                     (q0(k,i,j)*den0(k,i,j)-flux(k,i,j)/dxx(i,j)-x1*cf1*dd0(k,i,j-1))/den1(k,i,j) ))
              !print*,'i q= ', i,qn(i),line,q0(i),q0(i-1),q0(i+1)
          
	      !- for debug purposes only
              if(debug==on) QN(k,i,j)=(Q0(k,i,j)*DEN0(k,i,j)-FLUX(k,i,J)/DXX(I,J)-X1*CF1*DD0(k,I,j-1))/DEN1(k,i,j)
              !		This statement IGNORES monotonic limitations.
              !			   you should uncomment this line and run your
              !			   advection calculation wth constant initial
              !			   mixing ratios everywhere. If you have properly
              !			   implemented this subroutine constant mixing
	      !			   ratios should be maintiained
              ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
              !  print*,'i q= ', i,qn(i),line,q0(i),q0(i-1),q0(i+1)
              !  print*,FLUX(I),x1,cf1!,dd0(i-1),den1(i)
              !  print*,q0(i)-FLUX(I)-x1*cf1!,dd0(i-1),den1(i)
	      !      if(qn(i) < 99.999) pause

	      flux(k,i,j-1)=dxx(i,j)*(qn(k,i,j)*den1(k,i,j) - q0(k,i,j)*den0(k,i,j)) + flux(k,i,j)!Eq-9b
        END IF                                                  
     END DO
  END DO
  !
  ! If periodic boundary conditions are assumed, it is necessary
  !   to recalculate the updated mixing ratio at cell IDIM if there
  !   is inflow to that cell from the boundary between IDIM and 1
  !   Here these statements are commented out, but should be uncommented
  !   if this subroutine is needed for periodic boundary conditions,
  !   and then one of the calling arguements to the subroutine is IPERIOD
  !   which is set to "1" if you assume period boundary conditions
  !      IF(IPERIOD==1) THEN
  !      IF(U(1).LT.ZR0.AND.U(IDIM).LT.ZR0)
  !xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
  !     &  QN(IDIM)=(Q0(IDIM)*DEN0(IDIM)-FLUX(0)/DXX(IDIM)+FLUX(IDIM-1)/
  !     &                DXX(IDIM))/DEN1(IDIM)
  !      END IF
  !
ENDDO !- big loop x-z


END SUBROUTINE Advec3d_Y

!---------------------------------------------------------------------------
SUBROUTINE Advec3d_Z(m1,m2,m3, ia,iz,ja,jz,q0,u,den0,den1,dt,dxx,dd0,qn)
!-------------------------
  ! This subroutine calculates change in mixing ratio (Q0) during time
  !  step DT due to advection along a grid IDIM in length. Mixing ratios
  !  from host code (C) are loaded into Q0 array, which is updated to QN.
  !  Velocities (U) and fluxes (FLUX) are specified at cell FACES, having
  !  dimensions 0:IDIM. U, Q0, QN, DXX and FLUX indices defined here:
  !  Densities at beg, end time (DEN0, DEN1) defined in HOST CODE  
  !
  ! I grid->   |  1  |  2  |  I-1  |   I  |..   ..|  IDIM  | <- host grid
  ! U-array-> u(0)  u(1)  u(2)   u(i-1)  u(i)           u(IDIM)
  ! C-array->  | C(1)| C(2)| C(I-1)| C(I) |..   ..| C(IDIM)| mixing ratio
  ! DXX-arry-> | Dx1 | Dx2 | DxI-1 | DxI  |..   ..| DxIDIM |
  ! Density->  | Dd1 | Dd2 | DdI-1 | DdI  |..   ..| DdIDIM |
  !                 Q0 defined along 0 - IDIM+1 cells:
  !    |       | QN  | QN  |  QN   |  QN  |       |   QN   |        |
  !    |   Q0--|-Q0--|-Q0--|--Q0 --|--Q0--|..   ..|-- Q0 --|--Q0    |
  !    |    0  | 1   |  2  | I-1   |  I   |       |  IDIM  | IDIM+1 |
  !   lower BC |             <---   Q0 grid   --->         | upper BC
  !           Boundary conditions are stored in Q0 cells 0 & IDIM+1 
  !
  !  Input to this subroutine, provided in common /sub/, and the calling
  !  arguments to this subroutine:
  !     IDIM - #of grid cells being updated
  !     Q0(0:IDIM+1)- Initial mixing ratio along 1-D array, with two
  !                 additional boundary value mixing ratios padded into the
  !                 0th and IDIM+1 cell locations
  !     U(0:IDIM)- velocities BETWEEN grid cells (at the "higher-I" edges of
  !                each grid cell in the array, units consistent with DX, DT
  !     DEN0(IDIM)- Initial fluid density, which needs to be updated during
  !                 multi-dimensional calculations, as noted in Calling code
  !     DEN1(IDIM)- Updated fluid density, which needs to be updated during
  !                 multi-dimensional calculations, as noted in calling code
  !     DT-         time step- units consistent with U
  !     DXX(IDIM)-  Grid cell length along advection direction, Units
  !                   consistent with DT and U
  !     DD0(0:IDIM)- Initial fluid density flowing BETWEEN each grid cell
  !                  (remains constant for all dimensions at the initial
  !                  fluid density of the 1st dimension of a 2-3 D calculation
  !               one can use UPSTREAM density here (DD0(I)= RHO0(I) if u>0
  !               or DD0(I)= RHO0(I+1) if u<0) where RHO0 is the initial
  !               fluid density at the beginning of the 1st dimensional
  !               advection step of a 2 or 3 D advection calculation done one
  !               step at a time 
  !
  !  Output of this subroutine is an updated mixing ratio array QN(IDIM)
  !
  IMPLICIT none
  
  INTEGER,INTENT(IN)                   :: m1,m2,m3, ia,iz,ja,jz
  INTEGER                  :: idime
  REAL   ,INTENT(in),DIMENSION(m1,m2,m3) :: q0
  REAL   ,INTENT(in),DIMENSION(m1,m2,m3)   :: u
  REAL   ,INTENT(in),DIMENSION(m1,m2,m3)     :: den0
  REAL   ,INTENT(in),DIMENSION(m1,m2,m3)     :: den1
  REAL   ,INTENT(in)                   :: dt
  REAL   ,INTENT(in),DIMENSION(m1)     :: dxx
  REAL   ,INTENT(in),DIMENSION(m1,m2,m3)   :: dd0
  REAL   ,INTENT(OUT),DIMENSION(m1,m2,m3)    :: qn

  REAL,DIMENSION(m1)                 :: flux
  REAL,DIMENSION(m1)                   :: vcmax
  REAL,DIMENSION(m1)                   :: vcmin
  LOGICAL,DIMENSION(m1)            :: imxmn
  REAL,PARAMETER                       :: zr0=0.0
  
  INTEGER :: i,j,k
  REAL    :: cf,cf1,ck1,ck2,x1,x1n

  integer, parameter :: debug =0
  integer, parameter :: on = 1
  !

 idime = m1 ! z dir

  qn = q0

 !- big loop y-x
 DO j=ja,jz ; DO i=ia,iz
 
  
  imxmn=.false.
  ! imxmn(idime+1)=.false.
  ! Identify local max and min, specify mixing ratio limits at new time
  DO  k=2,idime-1 !
  
     imxmn(k)=q0(k,i,j)>=max(q0(k-1,i,j),q0(k+1,i,j)) .or. & !=true if local
              q0(k,i,j)<=min(q0(k-1,i,j),q0(k+1,i,j))	    !	    extrema
     ck1= q0(k,i,j)
     ck2= q0(k,i,j)
     if(u(k  ,i,j)< zr0) ck1= q0(k+1,i,j)
     if(u(k-1,i,j)>=zr0) ck2= q0(k-1,i,j)
     vcmax(k)= max( q0(k,i,j), ck1, ck2 )                      ! Eq-7
     vcmin(k)= min( q0(k,i,j), ck1, ck2 )                      ! Eq-7
  ENDDO
  !  VCMAX and VCMIN are the absolute physical limits to the
  !     mixing ratio at t+dt. If these limits are ever violated,
  !     non-monotonic (oscillatory) behavior in solution results
  !
  ! Update mixing ratios and limit Fluxes going UP where u>0
  !  First assume upstream flux at edge of domain
  IF(u(1,i,j)>=zr0) flux(1)= q0(1,i,j)*u(1,i,j)*dt*dd0(1,i,j)
  
  DO k=2,idime-1 ! ja,jz
  
  !print*,'part 1'
     IF(u(k,i  ,j)<zr0) CYCLE
     
     IF(u(k-1,i,j)<zr0) THEN
        flux(k)= q0(k,i,j)*u(k,i,j)*dt*dd0(k,i,j)    !  outflow-only cell
     
     ELSE                              !      use upstream
        x1= dt*u(k,i,j)/dxx(k)               ! Courant number
        x1n= (1.-x1)*(q0(k+1,i,j)-q0(k-1,i,j))/4.
        !
        ! First, estimate mixing ratio in outflowing fluid (Cf)
        cf= q0(k,i,j) + x1n                                       !Eq-4a
        !
        !   Check to see if there is a peak (min) upwind and/or
        !    downwind of cell face
        IF(imxmn(k-1)) cf= q0(k,i,j) +MAX(1.5,1.2  +.6 *x1)*x1n   !Eq-10b
        IF(imxmn(k+1)) cf= q0(k,i,j) +       (1.75 -.45*x1)*x1n   !Eq-10a
        !        CF= Q0(k,i,j) + 5.*X1N   ! uncomment this line for "full sharp"
        !
        !   Limit Cf to be between mixing ratio on either side of edge
        !      where flux is being calculated
        cf1= MIN( MAX( cf, MIN(q0(k,i,j),q0(k+1,i,j))  ), MAX(q0(k,i,j),q0(k+1,i,j)) )
        !- for debug purposes only
	if(debug==on)        CF1= CF     ! This statement IGNORES monotonic limitations.
        !                       you should uncomment this line and run your
        !                       advection calculation wth constant initial
        !                       mixing ratios everywhere. If you have properly
        !                       implemented this subroutine constant mixing
        !                       ratios should be maintained
        ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
        
	!
        !   Calculate mixing ratio at new time, but limit to physically
        !    reasonable values
        qn(k,i,j)= MAX(  vcmin(k),   MIN(   vcmax(k),          &   !eq-3&8
        (q0(k,i,j)*den0(k,i,j)-x1*cf1*dd0(k,i,j)+flux(k-1)/dxx(k))/den1(k,i,j) ))

	! for debug purposes only
        if(debug==on)  QN(k,i,j)=(Q0(k,i,j)*DEN0(k,i,j)-X1*CF1*DD0(k,i,j)+FLUX(k-1)/dxx(k))/DEN1(k,i,j)
        !            This statement IGNORES monotonic limitations.
        !                       you should uncomment this line and run your
        !                       advection calculation wth constant initial
        !                       mixing ratios everywhere. If you have properly
	!                       implemented this subroutine constant mixing
        !                       ratios should be maintiained
        ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
       !   print*,'i q= ', i,qn(i),line,q0(i),q0(i-1),q0(i+1)
	!       if(qn(i) < 99.999) pause

	!
        !   Re-calculate OUTFLOWING flux before moving on to next cell
        !    Flux = CF1*X1*DD0 but it must be adjusted if a monotonic limit
        !    is encountered.
        flux(k)= dxx(k)*(q0(k,i,j)*den0(k,i,j) - qn(k,i,j)*den1(k,i,j)) + flux(k-1)
     END IF                                                  !Eq-9a
  END DO
  !
  ! If periodic boundary conditions are assumed, it is necessary
  !   to recalculate the updated mixing ratio at cell 1 if there
  !   is inflow to that cell from the boundary between IDIM and 1
  !   Here these statements are commented out, but should be uncommented
  !   if this subroutine is needed for periodic boundary conditions,
  !   and then one of the calling arguements to the subroutine is IPERIOD
  !   which is set to "1" if you assume period boundary conditions
  !      IF(IPERIOD==1) THEN
  !        IF(U(IDIM-1)>=ZR0.AND.U(IDIM)>=ZR0)
  !     &  QN(1)=(Q0(1)*DEN0(1)-FLUX(1)/DXX(1)+FLUX(IDIM)/DXX(1))/DEN1(1)
  !      END IF 
  !
  ! Update mixing ratios and limit Fluxes going DOWN where u<0
  !  The logic of this loop through the grid line is identical
  !  to the "DO 10" Loop above, only you start at the highest I
  !  edge and work backwards to I=1
  !
  
  
  IF(u(idime-1,i,j)<zr0) flux(idime-1)=q0(idime,i,j)*u(idime-1,k,j)*dt*dd0(idime-1,i,j)
 
  DO k=idime-1,2,-1 !jz,ja,-1
     
       IF(u(k-1,i,j)>=zr0) THEN           ! Inflow-only cell
        
	IF(u(k,i,j)<zr0) qn(k,i,j)=  MAX(  vcmin(k),   MIN(   vcmax(k),&
           (q0(k,i,j)*den0(k,i,j)-flux(k)/dxx(k) + flux(k-1)/dxx(k))/den1(k,i,j) ))
        ELSE
           x1=  dt*ABS(u(k-1,i,j))/dxx(k)     ! Courant number
           x1n= (1.-x1)*(q0(k-1,i,j)-q0(k+1,i,j))/4.
           cf= q0(k,i,j) + x1n                                       !Eq-4b
           IF(imxmn(k+1)) cf= q0(k,i,j) +MAX(1.5,1.2  +.6 *x1)*x1n   !Eq-10b
           IF(imxmn(k-1)) cf= q0(k,i,j) +   (1.75 -.45*x1)*x1n       !Eq-10a
           cf1= MIN( MAX( cf, MIN(q0(k,i,j),q0(k-1,i,j)) ), MAX(q0(k,i,j),q0(k-1,i,j)) )
           !- for debug purposes only
	   if(debug==on)  CF1= CF	! This statement IGNORES monotonic limitations.
           !			   you should uncomment this line and run your
           !			   advection calculation wth constant initial
           !			   mixing ratios everywhere. If you have properly
           !			   implemented this subroutine constant mixing
           !			   ratios should be maintained
           ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
              
	   IF(u(k,i,j)>=zr0) cf1= q0(k,i,j)     ! outflow-only cell upstream
           qn(k,i,j)= MAX(  vcmin(k),  MIN(   vcmax(k), 	  &   !Eq-3&8
               (q0(k,i,j)*den0(k,i,j)-flux(k)/dxx(k)-x1*cf1*dd0(k-1,i,j))/den1(k,i,j) ))
           !print*,'i q= ', i,qn(i),line,q0(i),q0(i-1),q0(i+1)
          
	   !- for debug purposes only
           if(debug==on) QN(k,i,j)=(Q0(k,i,j)*DEN0(k,i,j)-FLUX(k)/dxx(k)-X1*CF1*DD0(k-1,I,j))/DEN1(k,i,j)
           !		This statement IGNORES monotonic limitations.
           !			   you should uncomment this line and run your
           !			   advection calculation wth constant initial
           !			   mixing ratios everywhere. If you have properly
           !			   implemented this subroutine constant mixing
	   !			   ratios should be maintiained
           ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
         !  print*,'i q= ', i,qn(i),line,q0(i),q0(i-1),q0(i+1)
         !  print*,FLUX(I),x1,cf1!,dd0(i-1),den1(i)
         !  print*,q0(i)-FLUX(I)-x1*cf1!,dd0(i-1),den1(i)
	 !      if(qn(i) < 99.999) pause

	   flux(k-1)=dxx(k)*(qn(k,i,j)*den1(k,i,j) - q0(k,i,j)*den0(k,i,j)) + flux(k)!Eq-9b
     END IF                                                  
  END DO
  !
  ! If periodic boundary conditions are assumed, it is necessary
  !   to recalculate the updated mixing ratio at cell IDIM if there
  !   is inflow to that cell from the boundary between IDIM and 1
  !   Here these statements are commented out, but should be uncommented
  !   if this subroutine is needed for periodic boundary conditions,
  !   and then one of the calling arguements to the subroutine is IPERIOD
  !   which is set to "1" if you assume period boundary conditions
  !      IF(IPERIOD==1) THEN
  !      IF(U(1).LT.ZR0.AND.U(IDIM).LT.ZR0)
  !xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
  !     &  QN(IDIM)=(Q0(IDIM)*DEN0(IDIM)-FLUX(0)/DXX(IDIM)+FLUX(IDIM-1)/
  !     &                DXX(IDIM))/DEN1(IDIM)
  !      END IF
  !
ENDDO;ENDDO !- big loop y-x


END SUBROUTINE Advec3d_Z

     
END MODULE monotonic_adv
